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INTRODUCTION 


Providing  accurate  medical  imaging  at  an  aid  station  or  remote  field  hospital  is  difficult. 
Portable  ultrasound  instrumentation  can  be  designed  for  these  applications,  but  the  expertise  required 
for  an  accurate  diagnosis  can  only  be  obtained  through  years  of  training.  Much  of  the  information 
in  an  ultrasound  examination  is  obtained  by  exploiting  the  real-time  nature  of  the  imaging  modality. 
A  successful  diagnosis  relies  on  the  skill  of  the  diagnostician  to  transform  mentally  dynamic  2-D 
images  into  the  complex  3-D  anatomy.  Locating  anatomical  landmarks  and  moving  the  scan  plane 
throughout  the  volume  of  interest  are  both  critical  components  of  this  process.  Without  extensive 
training,  it  would  be  difficult  for  a  medical  corpsman  to  perform  this  procedure. 

Researchers  have  started  exploring  tele-medicine  in  combination  with  three-dimensional  (3- 
D)  ultrasound  data  acquisition  as  a  solution  to  this  problem.  This  combination  could  potentially 
transfer  the  skill  required  to  scan  the  patient  and  make  the  diagnosis  from  a  medical  corpsman  to  an 
imaging  expert.  Unfortunately,  acquiring  a  3-D  ultrasound  data  set  is  difficult.  Modem  3-D 
ultrasound  systems  are  essentially  conventional  scanners  modified  to  collect  a  series  of  2-D  images. 
The  images  are  later  ‘stacked’  to  represent  the  3-D  anatomy.  Although  modem  scanners  are  designed 
to  collect  2-D  images  in  real-time  (20  2-D  images/s),  3-D  image  acquisition  is  slow.  Slow  image 
acquisition  introduces  the  problem  of  how  to  align  adjacent  2-D  images  collected  at  different  times. 
The  patient  and  imaging  probe  can  be  immobilized  to  reduce  movement  of  the  anatomy  between 
adjacent  images.  Cardiac  and  respiratory  gating  can  also  be  applied.  Even  in  a  carefully  controlled 
clinical  setting,  the  resulting  3-D  data  set  is  often  badly  distorted.  If  the  patient,  the  anatomy,  or  the 
transducer  moves  during  3-D  image  acquisition,  the  data  must  be  discarded.  Providing  the  care 
needed  to  obtain  ‘good’  3-D  data  in  the  clinic  is  troublesome,  on  a  battlefield  it  would  be  very 
difficult. 

We  are  building  an  ultrasound  imaging  system  that  would  avoid  these  problems.  The  two 
key  components  of  our  system  are  the  following:  1)  a  high  speed  scanner  that  can  collect  3-D  data 
40  to  80  times  faster  than  current  3-D  imaging  approaches,  and  2)  a  set  of  software  tools  for  rapid 
image  manipulation  and  analysis  at  a  remote  site.  Real-time  3-D  image  acquisition  eliminates  the 
need  for  patient  immobilization,  and  cardiac  and  respiratory  gating.  A  medical  corpsman  would 
simply  place  a  small  probe  on  the  patient  and  position  the  sample  volume  by  viewing  a  real-time 
two-dimensional  image  of  the  anatomy.  Once  the  probe  is  correctly  placed,  a  three-dimensional  data 
set  would  be  recorded  in  real-time  (0.05  s  for  each  3-D  data  set).  Following  data  acquisition,  the 
images  would  be  transmitted  to  a  central  hospital  for  post-processing  and  analysis.  An  expert 
clinician  could  ‘re-scan’  the  patient  by  looking  at  multiple  2-D  planes  through  the  data  sets,  or 
examine  a  computer  reconstruction  of  the  three-dimensional  anatomy.  The  3-D  data  set  could  also 
be  analyzed  quantitatively  to  calculate  dynamic  and/or  anatomic  changes  in  the  anatomy.  The  entire 
imaging  procedure,  including  subject  preparation,  would  take  only  a  few  minutes. 
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BODY 


The  Statement  of  Work  and  a  description  of  progress  in  each  area  is  given  below: 

Phase  I :  Month  1-6 
Scanner/Array  Development: 

•  System  specification  -  Perform  a  computer  simulation  to  evaluate  the  radiation  pattern  and 
imaging  performance  of  different  array  designs  and  beamforming  strategies.  Based  on  these 
results  derive  specifications  for  each  component  of  the  system. 

COMPLETED  -  This  work  is  described  in  the  publication:  “Theoretical  Assessment  of  a  Synthetic 
Aperture  Beamformer  for  Real-Time  3-D  Imaging”,  C.R.  Hazard  and  G.R.  Lockwood,  IEEE 
Transactions  on  Ultrasonics,  Ferroelectrics  and  Frequency  Control,  Vol.  4,  pp.  972-979, 1999.  (see 
Appendix  1) 

•  DSP  boards  -  Place  order  with  Pentek  Inc.  for  multi  processor  DSP  +  A/D  boards. 
COMPLETED  -  The  order  was  placed  03/25/99 

•  Transmit  electronics  -  Preliminary  design  study  for  high  voltage  transmit  circuit. 

NOT  COMPLETED  -  The  specifications  for  the  transmit  circuits  will  be  determined  by  the 
transducer  technology  (multi-layer  vs.  single  layer).  A  decision  on  the  transducer  technology  is 
expected  in  January  2000. 

•  Receive  electronics  -  Preliminary  design  study  for  front  end  electronics  (pre-amplifier,  time- 
gain-compensation  circuit  and  anti-aliasing  filter). 

COMPLETED  -  This  work  is  described  in  a  Masters  Thesis  by  Kurtis  Swinehart,  “Signal 
Conditioning  for  Real-Time  Three-Dimensional  Echocardiography”  Chapters  3-4,  Case  Western 
Reserve  University,  August,  1999.  (see  Appendix  2) 

The  following  tasks  were  assigned  to  Tetrad  Corporation  under  a  subcontract.  The  subcontract  was 
awarded  08/18/99.  A  report  describing  Tetrad’s  initial  progress  on  this  phase  of  the  work  was 
submitted  in  October,  99.  A  copy  of  this  report  is  included  in  Appendix  3. 

•  Transducer  development  -  Develop  methods  for  dicing  the  stacked  transducer  and  making 
interconnections. 

•  Verification  of  acoustic  output  capability-  Build  a  transducer  array  (several  elements)  and 
test  output  power. 

•  Phased  array  performance  optimization  -  Use  finite  element  and  1-D  transmission  line 
models  in  conjunction  with  partial  arrays  to  optimize  the  array  design  (matching  layers. 
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backing,  cut  geometry,  kerf  filler,  lens). 

•  Mechanical  probe  design  -  Preliminary  design  of  mechanical  probe,  identify  risk  areas  in 
mechanical  sector  scanner  and  perform  analysis  experiments  to  assure  the  required 
performance  can  be  met.  This  will  include  power  and  force  requirements,  rear  seal 
termination  method,  acoustic  window  and  fluid. 

•  Array  interconnection  -  Choose  a  method  for  connecting  the  array  elements  to  the  cable  that 
can  be  easily  packaged  and  will  tolerate  a  large  number  of  oscillations  at  10  Hz.  Test  and 
verify  approached  with  bench  testing. 

User  Interface  &  Image  Processing  Software: 

•  Develop  general  GUI  with  2D  and  3D  interactive  displays. 

COMPLETED  -  A  GUI  for  the  display  of  ultrasound  data  was  created  using  Tcl/Tk,  Togl,  C/C++ 
and  OpenGL  computer  languages  and  toolkits.  The  source  of  the  3D  ultrasound  images  has  been 
Volumetries,  Inc.’s  3D  ultrasound  scanner  located  in  the  Department  of  Cardiology  at  the  Cleveland 
Clinic  Foundation.  The  choice  of  the  programming  languages  and  toolkits  makes  the  GUI 
cross-platform.  New  controls  will  be  added  to  this  GUI  to  incorporate  new  processing  features. 

•  Develop  viewing  tool  for  3D  ultrasound  images. 

COMPLETED  -  Multi  Planar  Reformatting  (MPR)  is  one  of  the  most  common  modes  of  viewing 
three  and  higher  dimensional  data.  The  MPR  for  viewing  of  volume  data  has  been  implemented. 

•  Apply  texture  mapping  to  allow  interactive  panning  through  image  volumes. 

COMPLETED  -The  above-stated  MPR  viewing  mode  requires  computationally  intensive  trilinear 
interpolation  of  native  data  samples.  Texture  mapping  is  a  hardware-accelerated  approach  to 
performing  fast  trilinear  interpolation.  The  MPR  viewing  mode  has  been  implemented  using  texture 
mapping. 

Phase  II:  Month  7  - 12 
Scanner/Array  Development: 

•  Beamformer  -  configure  DSP  boards  and  write  and  test  assembly  language  code  for  single 
channel. 

COMPLETED  -  The  configuration  for  the  DSP  network  has  been  carefully  planned  to  account  for 
the  processor  to  processor  communication  capabilities  as  well  as  the  board  to  board  communication 
facilities  of  the  Pentek  4290  DSP  boards.  This  configuration  allows  for  optimal  communication 
among  the  processors  in  the  network.  Each  channel  in  the  system  has  one  TMS320C6201  processor 
which  performs  the  required  linear  interpolation  and  a  second  processor  which  is  involved  in  the  summing 
and  channel  to  channel  communications.  This  arrangement  has  been  evaluated  using  a  two  channel 
system.  Assembly  language  code  has  been  written  for  a  single  channel.  The  linear  interpolation 
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algorithm  has  been  tested  and  optimized.  This  code  is  capable  of  performing  the  required  linear 
interpolations  for  both  the  16  bit  inphase  (I)  and  16  bit  quadrature  (Q)  components  of  the  signal  at 
a  rate  of  one  I  and  Q  every  5  CPU  cycles.  Additional  code  has  been  written  to  test  the 
communication  algorithms  for  inter-processor  and  inter-board  communications.  All  of  the  testing 
has  been  done  on  an  initial  beta  version  of  the  DSP  board  which  operates  at  a  slower  speed  than  the 
actual  board  which  will  be  used. 

•  Build  and  test  high  voltage  transmit  electronics  (single  channel),  design  transmit  delay 
circuit. 

NOT  COMPLETED  -  A  subcontract  to  Queen’s  University  to  complete  this  work  was  awarded 
1 1/10/99.  We  are  still  waiting  for  a  decision  on  the  transducer  technology  (expected  Jan.  2000)  to 
proceed  with  this  work. 

•  Build  and  test  front  end  electronics  (single  channel). 

COMPLETED  -  See  Appendix  2. 

The  following  tasks  were  assigned  to  Tetrad  Corporation.  A  short  report  describing  Tetrad’s  progress 
over  the  last  quarter  is  attached  in  Appendix  4. 

•  Array  Design  -  based  on  results  of  Phase  I  complete  a  full  array  design. 

•  Purchase  parts  and  tooling  for  array  fabrication  and  built  initial  array  prototypes. 

•  Generate  a  design  for  mechanical  sector,  test  drive  electronics. 

•  Perform  life  testing  on  cables. 

•  Initial  acoustic  testing  with  array  in  probe  head. 

•  Specify  interface  requirements  for  the  array  and  drive  circuitry. 

User  Interface  &  Image  Processing  Software: 

•  Implement  surface  rendering  of  image  volumes. 

IN  PROGRESS  -  A  prerequisite  for  surface  rendering  of  image  volumes  is  that  surfaces  be  identified 
within  the  volume.  The  surface  identification  process  of  image  segmentation  is  underway.  We  are 
investigating  two  approaches.  The  first  semi-automatic  approach  based  on  deformable  models  uses 
a  priori  shape  information  to  capture  the  general  shape  of  the  structure  to  be  segmented.  This  model 
is  then  deformed  based  on  the  local  image  features  to  fit  the  true  underlying  shape.  The  second 
approach  is  an  extension  of  a  multi-resolution  bayesian  approach  to  segment  2D  ultrasound  images 
into  3D.  We  expect  our  initial  results  by  March  2000. 

•  Implement  volume  rendering  of  image  volumes. 

IN  PROGRESS  -  Volume  rendering  of  image  volumes  has  been  performed  using  a  specialized 
OpenGL-based  toolkit  called  Volumizer  and  a  program  called  VolView.  Similar  to  surface 
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rendering,  the  quality  of  volume  rendering  can  be  significantly  improved  when  the  result  of  image 
segmentation  is  available.  Image  segmentation,  which  is  underway  as  stated  above,  will  further  aid 
volume  rendering. 

•  Integrate  the  appropriate  controls  into  user  interface. 

COMPLETED  -  The  application  GUI  has  been  developed  with  extensibility  as  one  of  its  design 
features.  The  use  of  Tk  as  a  GUI  builder  language  allows  us  to  add  new  controls  into  the  user 
interface  conveniently.  The  controls  for  image  segmentation,  image  registration  and  the  proposed 
toolboxes  will  be  added  when  these  tasks  are  fully  developed. 


KEY  RESEARCH  ACCOMPLISHMENTS 

•  Theoretical  assessment  of  the  proposed  synthetic  aperture  beamformer  (see  Appendix  1) 

•  Preliminary  design  of  front  end  electronics  (see  Appendix  2) 

•  Identified  the  major  technical  risks  for  the  proposed  transducer  technology  (see  Appendix 
3) 

•  Evaluation  of  motion  artifacts  associated  with  the  proposed  synthetic  aperture  beamformer 
(see  Appendix  5) 

•  Interactive  multi-planar  reformatting 

•  Demonstration  of  registration  of  ultrasound  volumes  (see  Appendix  5) 


REPORTABLE  OUTCOMES 

Manuscripts,  abstracts,  presentations: 

i)  C.R.  Hazard  and  G.R.  Lockwood,  “Theoretical  assessment  of  a  synthetic  aperture  beamformer  for 
real-time  3-D  imaging,”  lEEETrans.  Ultrason.,  Ferroelect.,  Freq.  Contr.,  vol.  46,  pp.  972-979, 1999 

ii)  C.R.  Hazard  and  G.R.  Lockwood,  “Effects  of  motion  on  a  synthetic  aperture  beamformer  for  real¬ 
time  3D  ultrasound,”  Proc.  1999  IEEE  Ultrason.  Symp.,  Lake  Tahoe,  Nevada,  In  Press. 

iii)  K.A.  Swinehart,  “Signal  conditioning  for  real-time  three-dimensional  echocardiography”, 
Masters  Thesis,  Case  Western  Reserve  University,  1999. 

iv)  Vladimir  Zagrodsky,  Raj  Shekhar,  and  J.  Fredrick  Comhill,  “Mutual  information  based 
registration  of  cardiac  ultrasound  volumes,”  Proceedings  of SP1E  Medical  Imaging  2000,  San  Diego, 
California,  In  Press. 
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Degrees  obtained: 

Kurtis  Andrew  Swinehart,  Masters  of  Science,  Case  Western  Reserve  University,  August  1999. 


CONCLUSIONS 

We  have  developed  a  method  for  high  speed  ultrasound  imaging.  The  method  permits  the 
acquisition  of  40  to  80  images  in  the  time  normally  required  to  collect  a  single  image.  The  increased 
acquisition  speed  will  be  used  to  collect  a  3-D  data  set  in  real-time.  The  objective  of  this  proposal 
is  to  build  a  prototype  scanner  and  develop  the  software  tool  required  to  analyze  and  interpret  the 
3-D  data. 

A  major  challenge  associated  with  the  development  of  this  scanner  is  how  to  design  an 
ultrasound  beamformer  that  can  process  the  ultrasound  signals  at  the  required  frame  rates.  During 
the  first  year  of  this  project  we  proposed  a  suitable  beamformer  architecture  and  theoretically 
evaluated  the  performance  of  our  design  using  computer  simulations.  The  simulations  proved  the 
validity  of  our  approach.  They  also  provided  a  tool  for  deriving  accurate  specifications  for  each 
component  in  the  system.  Based  on  this  information  we  have  developed  a  detailed  design  for  the 
beamformer  and  have  started  the  process  of  assembling  the  electronic  hardware. 

Next  to  the  ultrasound  beamformer,  the  ultrasound  transducer  array  is  the  most  important  and 
difficult  component  to  design.  Based  on  our  specifications  for  the  array.  Tetrad  Corp.  has  started 
investigating  suitable  transducer  technologies  and  array  designs.  Preliminary  theoretical  work  and 
experimental  work  by  Tetrad  Corp.  have  suggested  that  a  new  multi-layer  transducer  design  should 
meet  or  exceed  our  specifications.  Unfortunately,  fabricating  a  multi-layer  device  is  significantly 
more  difficult,  and  issues  related  fabrication  and  reliability  are  still  being  investigated.  We  expect 
to  make  a  decision  about  the  transducer  technology  early  this  year  (2000). 

The  user  interface  and  image  analysis  software  development  is  underway  alongside  the 
scanner  development.  Currently,  the  source  of  3D  images  is  a  scanner  manufactured  by  Volumetries, 
Inc.  We  have  developed  software  tools  to  scan  convert  the  native  polar  data  for  image  analysis 
purposes.  The  development  of  interactive  MPR  viewing  mode  together  with  necessary  user  interface 
controls  is  complete.  We  have  also  performed  preliminary  volume  rendering  of  image  data. 

Ultrasound  image  segmentation  and  registration  algorithms  are  being  developed.  These 
algorithms  will  be  core  techniques  for  the  specialized  toolboxes  proposed  for  the  second  year  of  the 
project. 


Real-time  3-D  imaging  combined  with  tele-medicine  and  a  set  of  image  analysis  tools  with 
enable  ultrasound  imaging  for  forward  echelon  combat  casualty  care.  This  technology  will  be  equally 
effective  in  civilian  emergency  care. 
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Abstract — A  real-time  3-D  imaging  system  requires  the 
development  of  a  beamformer  that  can  generate  many 
beams  simultaneously.  In  this  paper,  we  discuss  and  eval¬ 
uate  a  suitable  synthetic  aperture  beamformer.  The  pro¬ 
posed  beamformer  is  based  on  a  pipelined  network  of  high 
speed  digital  signal  processors  (DSP).  By  using  simple 
interpolation-based  beamforming,  only  a  few  calculations 
per  pixel  are  required  for  each  channel,  and  an  entire  2-D 
synthetic  aperture  image  can  be  formed  in  the  time  of  one 
transmit  event.  The  performance  of  this  beamformer  was 
explored  using  a  computer  simulation  of  the  radiation  pat¬ 
tern.  The  simulations  were  done  for  a  full  64-element  array 
and  a  sparse  array  with  the  same  receive  aperture  but  only 
five  transmit  elements.  We  assessed  the  effects  of  changing 
the  sampling  rate  and  amplitude  quantization  by  compar¬ 
ing  the  relative  levels  of  secondary  lobes  in  the  radiation 
patterns.  The  results  show  that  the  proposed  beamformer 
produces  a  radiation  pattern  equivalent  to  a  conventional 
beamformer  using  baseband  demodulation,  provided  that 
the  sampling  rate  is  approximately  10  times  the  center  fre¬ 
quency  of  the  transducer  (34%  bandwidth  pulse).  The  sim¬ 
ulations  also  show  that  the  sparse  array  is  not  significantly 
more  sensitive  to  delay  or  amplitude  quantization  than  the 
full  array. 


I.  Introduction 

Real-time  3-D  ultrasound  imaging  has  the  potential  to 
provide  unique  clinical  information  for  the  diagnosis 
and  treatment  of  disease.  We  recently  described  a  method 
for  real-time  3-D  imaging  using  a  sparse  linear  phased  ar¬ 
ray  and  synthetic  aperture  beamforming  [1],  [2].  A  real¬ 
time  3-D  data  set  is  collected  by  mechanically  rocking 
the  transducer  array  in  the  elevation  direction  and  rapidly 
acquiring  2-D  images.  Computer  simulations  showed  that 
frame  rates  as  high  as  one  thousand  two-dimensional  im¬ 
ages  per  second  can  be  achieved  using  synthetic  aperture 
beamforming  without  a  reduction  in  the  signal-to-noise  ra¬ 
tio  and  with  little  degradation  of  image  contrast  and  reso¬ 
lution.  However,  developing  this  real-time  3-D  system  will 
require  a  beamformer  capable  of  generating  many  beams 
simultaneously.  Most  modern  beamformers  use  a  dedicated 
beamformer  channel  for  each  active  receive  element.  This 
allows  parallel  processing  of  the  received  signals  and  real¬ 
time  beamforming.  Because  only  one  beam  is  generated 
at  a  time,  the  same  beamformer  hardware  can  be  used 
for  subsequent  lines  in  the  image  by  simply  varying  the 
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delay  and  apodization  pattern.  This  approach  is  not  prac¬ 
tical  for  a  real-time  3-D  imaging  system.  If  a  conventional 
beamformer  were  simply  duplicated  to  generate  the  beams 
required  for  real-time  3-D  imaging,  the  resulting  system 
would  be  massive.  A  beamformer  for  a  64-element  array 
generating  80  receive  beams  simultaneously  would  require 
5120  beamformer  channels.  Some  reduction  in  hardware 
can  be  achieved  by  separating  the  expensive  focusing  de¬ 
lays  and  the  fixed  steering  delays.  This  allows  the  same 
focusing  delays  to  be  shared  by  adjacent  beams  provided 
the  angular  spacing  is  not  too  great.  This  innovative  ap¬ 
proach  was  first  introduced  by  von  Ramm  et  al  [3]  and 
Shattuck  et  al.  [4]  and  was  used  for  a  real-time  3-D  imag¬ 
ing  system  based  on  a  2-D  transducer  array.  Unfortunately, 
this  technique  is  not  appropriate  for  the  proposed  syn¬ 
thetic  aperture  system  because  the  beams  must  be  formed 
simultaneously  throughout  the  entire  image. 

The  need  to  generate  multiple  receive  beams  simulta¬ 
neously  can  be  avoided  by  storing  the  signals  from  each 
element  in  random  access  memory.  Once  the  signals  are 
stored,  an  image  can  be  formed  by  adding  points  from 
memory  locations  selected  to  provide  the  samples  closest 
to  the  desired  delayed  values.  This  nearest  neighbor  ap¬ 
proach  to  beamforming  has  three  advantages.  First,  the 
complexity  of  the  beamformer  is  reduced  because  only 
those  points  that  contribute  directly  to  each  pixel  are  used. 
Second,  the  pixels  can  be  generated  on  a  rectangular  grid 
and  displayed  without  a  separate  scan  conversion.  Finally, 
and  most  importantly,  storing  the  digitized  data  in  a  ran¬ 
dom  access  memory  allows  serial  processing  of  the  pixels. 
Serial  processing  avoids  the  hardware  duplication  associ¬ 
ated  with  parallel  beam  formation.  However,  for  real-time 
3-D  imaging  the  processing  speed  must  be  increased  in 
proportion  to  the  number  of  beams.  Corl  and  Kino  [6]  and 
others  [5],  [7],  [8]  developed  a  nearest  neighbor  beamformer 
for  a  2-D  synthetic  aperture  imaging  system  in  1978.  Be¬ 
cause  of  hardware  limitations,  the  signals  were  sampled  at 
close  to  the  Nyquist  frequency,  and  only  a  single  receive 
signal  was  recorded  for  each  transmit  event.  The  resulting 
delay  accuracy  of  the  system  was  too  coarse  to  provide 
adequate  suppression  of  secondary  lobes  in  the  radiation 
pattern,  and  the  images  were  not  suitable  for  medical  ap¬ 
plications. 

Modern  digital  beamformers  avoid  the  problem  of  ex¬ 
cessive  sampling  by  using  a  combination  of  coarse  and  fine 
delays.  Coarse  delays  are  determined  by  the  sampling  pe¬ 
riod,  and  fine  delays  are  generated  within  the  beamformer. 
Many  different  methods  have  been  suggested  for  generat¬ 
ing  fine  delays  based  on  up-sampling  or  interpolating  the 
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data,  and  a  variety  of  clever  hardware  designs  have  been 
proposed  [9]— [13] .  Although  up-sampling  or  interpolation 
reduces  the  speed  required  for  the  analog- to- digital  con¬ 
verters  (A/D),  the  number  of  points  that  must  be  gen¬ 
erated  and  processed  by  the  beamformer  is  still  large.  A 
common  alternative  to  up-sampling  is  to  heterodyne  the 
signals  to  an  intermediate  frequency  or  baseband  before 
beamforming  [14]— [23].  Simple  delay  and  add  processing 
at  data  rates  much  closer  to  the  Nyquist  frequency  can 
then  be  used,  provided  the  phase  of  each  signal  is  adjusted 
appropriately.  Heterodyning  reduces  the  data  rate,  but 
additional  costs  associated  with  mixing  and  filtering  the 
signals  may  negate  this  advantage.  Newer  beamforming 
techniques  based  on  oversampling  delta-sigma  A/D  offer 
promise  for  the  future  [24]— [26],  but  current  digital  beam- 
formers  are  still  large  and  expensive.  Excellent  reviews  of 
digital  beamforming  techniques  are  given  by  Thomenius 
[27]  and  Mucci  [28]. 

In  this  paper,  we  investigate  a  synthetic  aperture  beam- 
former  architecture  suitable  for  real-time  3-D  imaging.  The 
beamformer  is  similar  to  the  synthetic  aperture  beam- 
former  proposed  by  Corl  and  Kino  with  the  addition  of 
a  fine  delay  correction  and  use  of  receive  signals  from  all 
elements  for  each  transmit  event.  We  show  that  these  mod¬ 
ifications  allow  high  quality  beam  formation  with  little  in¬ 
crease  in  the  sampling  rate  and  only  a  moderate  increase  in 
the  complexity  of  the  beamformer.  The  performance  of  the 
proposed  beamformer  as  a  function  of  the  sampling  rate 
and  amplitude  quantization  was  evaluated  using  a  com¬ 
puter  simulation  of  the  radiation  pattern.  The  results  are 
compared  with  those  for  a  nearest  neighbor  beamformer 
and  a  modern  beamformer  using  baseband  interpolation. 


II.  Synthetic  Aperture  Beamformer  Design 

The  beamformer  proposed  by  Corl  and  Kino  [6]  for  2- 
D  synthetic  aperture  imaging  has  many  features  that  are 
attractive  for  a  real-time  3-D  imaging  system.  However, 
the  technology  available  in  1978  limited  the  sampling  fre¬ 
quency,  did  not  allow  multiple  receive  channels,  and  pre¬ 
vented  real-time  interpolation  of  the  synthetic  aperture 
data.  For  a  medical  ultrasound  beamformer,  interpolation 
is  essential  to  provide  the  fine  delays  required  for  adequate 
suppression  of  secondary  lobes  and  to  avoid  excessive  over- 
sampling.  Unfortunately,  the  cost  of  interpolating  the  data 
for  a  real-time  3-D  system  is  large.  For  example,  3.2  billion 
interpolations  per  second  would  be  required  for  real-time 

0.2  ms)  generation  of  the  inphase  and  quadrature  com¬ 
ponents  associated  with  a  5000  pixel  image  made  using 
a  64-element  array.  Obviously,  the  number  of  calculations 
required  for  each  interpolated  point  must  be  minimized. 
Another  problem  associated  with  real-time  interpolation 
of  a  synthetic  aperture  data  set  is  that  the  samples  are 
not  processed  in  the  order  received.  Samples  may  be  used 
for  more  than  one  pixel  or  not  used  at  all.  A  synthetic 
aperture  beamformer  must  have  the  flexiblity  to  adjust 
the  order  in  which  samples  are  processed.  This  sorting  re¬ 


quirement  makes  it  difficult  and  inefficient  to  use  a  con¬ 
ventional  up-sampling  or  heterodyning  approach  for  fine 
delay  generation. 

DSP  offer  an  attractive  solution  to  these  problems.  DSP 
provide  on-chip  memory  storage  and  the  memory  manage¬ 
ment  needed  to  sort  the  data  for  each  pixel.  Modern  DSP 
also  provide  the  processing  speed  necessary  for  real-time 
interpolation.  Fig.  1  shows  a  block  diagram  of  the  proposed 
beamformer.  The  beamformer  has  two  parallel  chains  with 
a  pair  of  DSP  associated  with  each  array  element.  The 
memory  inside  each  DSP  is  arranged  into  three  separate 
blocks.  One  block  is  used  to  store  the  digitized  data,  and 
two  blocks  are  organized  as  stacks  that  contain  the  mem¬ 
ory  addresses  corresponding  to  the  coarse  delays  and  the 
coefficients  for  interpolation  and  apodization.  The  coarse 
delays  and  interpolation  and  apodization  coefficients  are 
calculated  ahead  of  time  and  loaded  into  memory  when 
the  system  is  turned  on.  Because  data  can  be  randomly 
accessed  from  memory,  each  DSP  can  work  on  a  different 
pixel,  and  the  order  of  calculations  can  be  arranged  for 
pipelined  addition. 

Fig.  2  illustrates  the  flow  of  data  and  the  timing  of  cal¬ 
culations  for  the  proposed  beamformer.  The  calculation 
of  each  pixel  begins  with  the  first  DSP  at  the  top  of  the 
pipeline  and  continues  through  the  chain  until  signals  from 
each  array  element  have  been  added.  For  each  calculation 
cycle,  each  DSP  must  calculate  an  interpolated  value  (; pyc , 
the  interpolated  value  for  pixel  p  at  channel  c)  and  add  it 
to  an  accumulating  sum  ( PSC ,  the  sum  for  pixel  p  at  chan¬ 
nel  c).  The  interpolated  value  is  calculated  by  first  fetching 
two  data  points  from  memory  that  are  closest  to  the  de¬ 
sired  delay  value.  A  fine  delay  correction  is  introduced  by 
linearly  interpolating  between  the  two  points.  Once  this 
calculation  is  complete,  the  interpolated  value  is  added  to 
the  accumulating  sum  that  was  passed  down  from  the  up¬ 
stream  channels.  The  updated  sum,  PSC ,  is  then  passed  to 
the  next  DSP.  If  V\  and  v 2  represent  the  sample  values 
surrounding  the  desired  delayed  value,  the  new  value  of 
the  accumulating  sum,  PSC ,  is  given  by: 

p£c  =  pyc  +  pSc-i  =  7((v2  -  vi)At  +  vl)  +  PSC 

(1) 

where  7  is  the  apodization  coefficient,  At  is  the  fraction  of 
the  sample  period  corresponding  to  the  desired  fine  delay, 
and  pSc_i  is  the  accumulating  sum  passed  down  from  the 
upstream  DSP.  Many  DSP  are  optimized  to  calculate  a 
multiply  and  add  operation  in  a  single  cycle.  Consequently, 
it  is  useful  to  rearrange  (1)  into  two  multiply  and  add 
operations: 

»^temp  =::  otv\  -{-  pSq—i  and  (2) 

PSc  =  pv2  +  Stemp  (3) 

where  a  =  7(1  —  At)  and  /?  =  7 (At).  Assuming  the  coeffi¬ 
cients  a  and  /?  are  calculated  ahead  of  time  and  stored  in 
memory,  only  two  multiply  and  add  operations  per  chan¬ 
nel  are  required  to  interpolate  and  apodize  the  inphase  or 
quadrature  component  of  a  pixel.  In  comparison,  a  simple 
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Fig.  1.  Proposed  synthetic  aperture  beamformer  architecture.  The  signals  received  by  each  element  are  digitized  and  stored  in  memory.  A 
network  of  DSP  is  used  to  delay  signals  using  linear  interpolation.  The  network  of  DSP  is  arranged  in  a  pipeline  architecture  to  facilitate 
summing  of  the  delayed  signals. 


20-tap  FIR  interpolation  filter  would  require  10  times  as 
many  calculations.  Another  advantage  of  this  approach  is 
that  different  values  of  a  and  /3  can  be  stored  in  each  DSP 
for  every  pixel  in  the  image.  This  provides  dynamic  con¬ 
trol  over  the  aperture  shape  and  delay  pattern  for  every 
transmit-receive  event  and  optimal  focusing  at  every  pixel. 

The  final  task  of  the  beamformer  is  to  detect  the  en¬ 
velope  of  the  signal.  Second-order  sampling  provides  an 
efficient  estimate  of  the  envelope  at  a  given  pixel  [17],  [29]. 
Parallel  processing  streams  are  used  to  calculate  two  in¬ 
terpolated  values  at  each  array  element.  One  value,  called 
the  inphase  component  (/),  is  calculated  for  the  desired  de¬ 
lay,  and  a  second  beamformed  value,  called  the  quadrature 
component  (Q),  is  calculated  with  an  additional  delay  of 


one-quarter  of  a  period.  Assuming  that  the  envelope  varies 
slowly,  the  envelope  can  be  approximated  as  the  square 
root  of  the  sum  of  the  squares: 

Envelope  =  i//2  +  Q2.  (4) 

A  final  DSP  at  the  end  of  the  pipeline  is  used  to  add  the 
I  and  Q  values  from  multiple  transmit  events  and  to  cal¬ 
culate  the  envelope  value  for  each  pixel.  Because  second- 
order  sampling  requires  that  two  interpolated  values  be 
calculated  for  each  pixel  for  each  receive  channel,  the  to¬ 
tal  cost  per  channel  for  calculating  a  pixel  is  four  multiply 
and  add  operations,  which  is  small  enough  to  permit  real¬ 
time  3-D  beamforming  using  existing  DSP  technology.  For 
example,  the  fixed  point  TMS320C6x  manufactured  by 
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Fig.  2.  Diagram  showing  the  flow  of  data  and  the  timing  for  the  calculation  of  a  pixel  using  the  proposed  linear  interpolation  beamformer. 
pyc  is  the  interpolated  contribution  to  pixel  p  for  channel  c,  and  PSC  is  the  accumulating  sum  for  pixel  p  at  channel  c. 


Texas  Instruments  (Dallas,  TX)  can  process  400  million 
multiply  and  add  operations  every  second. 

Linear  interpolation  provides  an  efficient  method  for 
generating  a  fine  delay,  provided  the  sampling  frequency 
is  adequate.  The  success  of  this  approach  will  ultimately 
depend  on  how  fast  the  data  must  be  sampled  for  sup¬ 
pression  of  secondary  lobes  (less  than  —60  dB).  Many  re¬ 
searchers  have  investigated  the  impact  of  amplitude  and 
delay  quantization  on  secondary  lobes  [8],  [14],  [15],  [22], 
[27],  [28],  [30]— [36] .  A  variety  of  system  parameters,  includ¬ 
ing  the  bandwidth  of  the  pulse;  the  number,  spacing,  and 
apodization  of  the  transmit  and  receive  elements;  and  the 
beamforming  algorithm,  have  all  been  shown  to  affect  the 
amplitude  of  secondary  lobes.  Generally,  there  is  an  inverse 
relationship  between  the  complexity  of  the  beamformer 
and  the  required  sampling  frequency.  Unfortunately,  ac¬ 
curately  predicting  the  amplitude  of  secondary  lobes  for 
a  given  system  is  difficult.  In  the  next  section,  we  use  a 
computer  simulation  to  investigate  the  performance  of  the 
proposed  beamformer  as  a  function  of  sampling  frequency 
and  amplitude  quantization. 

III.  Computer  Simulations 

The  effects  of  amplitude  quantization  and  discrete  time 
sampling  on  the  radiation  pattern  of  an  array  were  exam¬ 
ined  for  three  synthetic  aperture  beamformers:  the  nearest 
neighbor  approach  introduced  by  Corl  and  Kino  [5]-[8],  the 


proposed  linear  interpolation  beamformer,  and  a  baseband 
demodulation  beamformer  using  the  method  outlined  by 
O’Donnell  et  al.  [14].  The  nearest  neighbor  approach  is 
the  simplest  and  least  expensive  computationally,  but  a 
prohibitively  high  sampling  rate  is  needed  for  adequate 
delay  accuracy.  At  the  other  extreme,  baseband  demodu¬ 
lation  requires  a  modest  sampling  rate  (4  to  10  times  the 
bandwidth  of  the  pulse),  but  the  required  mixing,  filter¬ 
ing,  and  phase  rotation  operations  associated  with  shifting 
the  signal  to  baseband  are  computationally  expensive.  The 
proposed  linear  interpolation  beamformer  offers  a  compro¬ 
mise  between  the  nearest  neighbor  and  heterodyning  ap¬ 
proaches  in  the  sampling  rate  and  the  number  of  compu¬ 
tations  per  pixel. 

For  each  of  the  three  different  beamforming  approaches, 
two  different  transducer  arrays  were  simulated:  a  full  array 
using  64  transmit  and  receive  elements  (0.95  A/2  element 
spacing)  and  a  sparse  array  with  the  same  receive  aper¬ 
ture  (64  elements)  but  only  5  active  transmit  elements 
[transmit  spacing  11  (0.95) A/2].  The  center  frequency  of 
the  transducer  array  was  3.5  MHz,  and  a  raised  cosine 
apodization  was  used  for  both  arrays.  The  speed  of  a  syn¬ 
thetic  aperture  imaging  system  is  determined  ultimately 
by  the  number  of  transmit  elements.  By  using  only  a  few 
transmit  events  per  image,  the  frame  rates  necessary  for 
real-time  3-D  imaging  can  be  obtained.  The  signal  loss 
associated  with  fewer  transmit  pulses  can  be  reduced  by 
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increasing  the  power  delivered  to  each  transmit  element 
and  by  using  multiple  transmit  elements  for  each  trans¬ 
mit  pulse  [1],  [37].  Although  a  full  array  with  64  transmit 
elements  would  be  too  slow  to  acquire  a  3-D  synthetic 
aperture  data  set  in  real-time,  the  radiation  pattern  of  the 
full  array  provides  a  useful  benchmark  for  comparing  the 
performance  of  the  sparse  array. 

The  broadband,  two-way  radiation  pattern  for  each  ar¬ 
ray  was  evaluated  numerically  by  simulating  a  sector  scan 
of  a  point  target.  Attenuation  was  ignored  in  this  sim¬ 
ulation,  and  a  uniform  element  directivity  was  assumed. 
The  point  target  was  placed  in  front  of  the  array  at  a  dis¬ 
tance  equal  to  four  times  the  receive  aperture  width.  Sig¬ 
nals  received  from  the  point  target  were  calculated  for  ev¬ 
ery  transmitter- receiver  pair  and  then  discretely  sampled 
in  time  and  quantized  in  amplitude  to  represent  analog- 
to-digital  conversion  at  the  front  end  of  the  beamformer. 
Radial  scan  lines  were  formed  from  -45°  to  45°  in  steps  of 
1°.  The  amplitude  of  the  signal  at  each  point  in  the  sector 
scan  was  calculated  by  summing  an  appropriately  delayed 
and  weighted  pulse  over  each  transmit  and  receive  element 
and  then  approximating  the  envelope  using  second-order 
sampling.  The  pulse  was  a  Gaussian-modulated  sinusoid 
with  a  center  frequency  of  3.5  MHz  and  a  —6  dB  band¬ 
width  of  34%.  Beamforming  delays  were  introduced  using 
one  of  the  three  beamforming  approaches.  We  repeated  the 
simulations  using  different  sampling  frequencies  (from  20 
to  100  MHz)  and  varying  levels  of  amplitude  quantization 
(from  8  to  14  bits).  The  radiation  patterns  were  displayed 
by  plotting  the  maximum  value  of  the  demodulated  signal 
along  each  scan  line  as  a  function  of  the  scan  angle.  This 
ensures  that  the  worst  case  secondary  lobes  are  displayed. 

The  effects  of  amplitude  and  delay  quantization  were 
considered  separately.  In  section  A,  the  radiation  patterns 
were  calculated  with  varying  levels  of  amplitude  quantiza¬ 
tion,  but  no  delay  quantization.  In  section  B,  delay  quan¬ 
tization  is  introduced  with  no  amplitude  quantization. 

A.  Amplitude  Quantization 

Fig.  3  shows  the  radiation  patterns  for  the  sparse  ar¬ 
ray  with  different  levels  of  amplitude  quantization.  Perfect 
beamforming  (no  amplitude  quantization  or  sampling)  is 
also  included  as  a  baseline  for  comparison.  The  main  lobe 
of  the  pattern  is  identical  to  the  perfect  beamforming  case 
for  all  levels  of  quantization  down  to  —60  dB.  For  8-bit 
quantization,  the  secondary  lobes  are  around  —70  dB;  for 
10  bits,  they  are  approximately  —80  dB;  at  12  bits,  they 
are  suppressed  to  -90  dB.  At  14  bits,  the  radiation  pat¬ 
tern  is  similar  to  the  perfect  beamforming  case;  therefore, 
the  secondary  lobes  are  limited  by  the  number  of  array 
elements  and  not  the  amplitude  quantization. 

Fig.  4  shows  the  effect  of  the  different  levels  of  ampli¬ 
tude  quantization  on  the  full  array.  The  main  lobe  is  again 
identical  to  the  perfect  beamforming  case  down  to  —60  dB. 
By  increasing  the  number  of  transmit  elements  from  5  to 
64,  the  level  of  quantization-limited  secondary  lobes  is  im¬ 
proved  by  approximately  5  to  10  dB.  For  the  full  array,  the 


Azimuthal  Angle  (deg.) 

Fig.  3.  Normalized  radiation  patterns  for  the  sparse  array  (5  transmit 
elements,  64  receive  elements)  using  8-,  10-,  12- ,  and  14-bit  quanti¬ 
zation.  The  unquantized  radiation  pattern  is  also  shown  for  compar¬ 
ison. 


Fig.  4.  Normalized  radiation  patterns  for  the  full  array  (64  transmit 
and  receive  elements)  using  8-,  10-,  12-,  and  14-bit  quantization.  The 
unquantized  radiation  pattern  is  also  shown  for  comparison. 

radiation  pattern  with  14-bit  quantization  is  similar  to  the 
perfect  beamforming  pattern  down  to  —110  dB. 

B.  Sampling  Rate 

The  effects  of  sampling  rate  on  secondary  lobes  were 
explored  for  the  three  beamforming  strategies:  nearest 
neighbor,  linear  interpolation,  and  baseband  demodula¬ 
tion.  Each  interpolation  method  was  evaluated  at  sam¬ 
pling  rates  from  20  to  100  MHz.  Fig.  5  shows  the  result¬ 
ing  secondary  lobe  levels  for  the  sparse  array  using  nearest 
neighbor,  linear  interpolation,  and  baseband  demodulation 
beamforming.  Sampling  rates  of  20,  60,  and  100  MHz  are 
shown  in  Figs.  5(a)-(c),  respectively.  Figs.  6(a)-(cj  show 
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Fig.  5.  Sparse  array  radiation  patterns  for  nearest  neighbor,  linear  in¬ 
terpolation,  and  heterodyning  (baseband  demodulation)  beamform- 
ers  with  sampling  rates  of  (a)  20,  (b)  60,  and  (c)  100  MHz. 


the  corresponding  results  for  the  full  array  with  64  trans¬ 
mit  elements. 

Using  nearest  neighbor  beamforming,  the  secondary 
lobes  for  the  sparse  array  (Fig.  5)  range  from  about  —35  dB 
at  a  20  MHz  sampling  rate  to  approximately  —50  dB  at 
100  MHz.  The  results  for  the  full  array  (Fig.  6)  are  only 
slightly  better  with  secondary  lobe  levels  from  —35  to 
-60  dB  over  the  same  range  of  sampling  frequencies.  A 
prohibitive  sampling  rate  of  about  400  MHz  is  required  to 
achieve  secondary  lobe  levels  below  —60  dB  using  near¬ 
est  neighbor  beamforming.  At  a  sampling  rate  of  20  MHz 
[Fig.  5(a),  6(a)],  baseband  demodulation  and  linear  inter¬ 
polation  beamforming  produce  similar  radiation  patterns 
for  both  the  full  and  sparse  arrays.  At  higher  sampling 
rates,  linear  interpolation  beamforming  produces  lower 
secondary  lobes  than  baseband  demodulation  beamform¬ 
ing.  There  is  a  much  greater  improvement  in  the  radiation 
pattern  as  a  function  of  sampling  rate  for  the  linear  inter¬ 
polation  beamformer  than  for  the  baseband  demodulation 
beamformer.  At  sampling  rates  greater  than  20  MHz,  the 
secondary  lobes  are  below  -60  dB  for  linear  interpolation 
and  baseband  beamforming.  The  sparse  array  was  not  sig¬ 
nificantly  more  sensitive  to  changes  in  the  sampling  rate 
than  the  full  array.  This  was  true  for  both  linear  interpo¬ 
lation  and  baseband  demodulation. 

Fig.  7  shows  the  radiation  pattern  for  the  sparse  ar¬ 
ray  using  linear  interpolation  beamforming  with  a  sam¬ 
pling  rate  of  40  MHz  and  10-bit  amplitude  quantization. 
Combining  the  effects  of  amplitude  quantization  and  sam¬ 
pling  further  degrades  the  radiation  pattern,  but  secondary 
lobes  are  still  below  —60  dB  at  all  angles. 


IV.  Summary  and  Conclusions 

In  this  paper,  we  described  a  beamformer  suitable  for 
real-time  3-D  synthetic  aperture  imaging.  The  beamformer 
is  realized  using  a  network  of  digital  signal  processors. 
Modern  DSP  provide  the  speed  and  the  memory  man¬ 
agement  required  for  sorting  and  interpolating  the  data 
in  real  time.  By  controlling  the  order  in  which  pixels  are 
generated,  the  DSP  network  can  be  arranged  in  a  pipeline 
structure  to  allow  efficient  summing  of  the  contributions 
from  each  receive  channel.  A  unique  set  of  interpolation 
coefficients  is  used  to  calculate  the  interpolated  samples 
required  for  each  pixel  in  the  image.  By  precalculating 
the  coefficients  and  using  second-order  sampling  for  en¬ 
velope  detection,  only  four  multiply  and  add  operations 
per  channel  are  required  for  each  pixel,  independent  of  the 
sampling  rate.  This  strategy  provides  an  accurate  and  very 
efficient  estimate  of  the  interpolated  samples  without  the 
proliferation  of  samples  associated  with  up-sampling  or  the 
complexity  of  a  heterodyning  approach. 

The  success  of  the  proposed  linear  interpolation  beam- 
former  will  depend  on  the  sampling  frequency  and  am¬ 
plitude  quantization  levels  needed  to  suppress  secondary 
lobes  adequately.  A  real-time  3-D  synthetic  aperture  sys¬ 
tem  requires  the  use  of  a  sparse  array.  Because  fewer 
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Fig.  6.  Full  array  (64  transmit  and  receive  elements)  radiation  pat¬ 
terns  for  nearest  neighbor,  linear  interpolation,  and  heterodyning 
(baseband  demodulation)  beamformers  with  sampling  rates  of  (a)  20, 
(b)  60,  and  (c)  100  MHz. 
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Fig.  7.  The  normalized  radiation  pattern  for  the  sparse  array  (5  trans¬ 
mit  elements,  64  receive  elements)  digitizing  at  40  MHz  with  10  bits 
of  accuracy  and  using  linear  interpolation  beamforming. 

transmit-receive  events  are  used  in  beamforming,  it  is  par¬ 
ticularly  important  to  evaluate  the  sensitivity  of  the  sys¬ 
tem  to  amplitude  and  delay  quantization.  The  effects  of 
sampling  and  amplitude  quantization  were  examined  using 
computer  simulations  of  the  beamformer.  The  simulations 
assumed  a  Gaussian-shaped  pulse,  centered  at  3.5  MHz 
with  a  -6  dB  bandwidth  of  34%.  Similar  results  were  ob¬ 
tained  for  a  pulse  with  a  —6  dB  bandwidth  of  68%.  The 
results  show  that  a  10-bit  A/D  would  provide  sufficient  dy¬ 
namic  range  to  keep  secondary  lobes  below  —  70  dB  for  the 
sparse  array  (5  transmit  elements,  64  receive  elements).  At 
finer  levels  of  amplitude  quantization,  secondary  lobes  are 
limited  by  the  number  of  array  elements  and  not  by  the 
amplitude  quantization.  The  secondary  lobes  in  the  radi¬ 
ation  pattern  of  the  linear  interpolation  beamformer  were 
found  to  be  suppressed  approximately  20  dB  compared 
with  the  nearest  neighbor  beamformer.  At  a  sampling  rate 
of  20  MHz,  baseband  demodulation  and  linear  interpola¬ 
tion  have  similar  performance.  At  higher  sampling  rates, 
linear  interpolation  outperforms  baseband  demodulation. 
The  sparse  array  was  not  found  to  be  significantly  more 
sensitive  to  amplitude  and  delay  quantization  than  the  full 
array  (64  transmit  and  receive  elements).  The  difference 
between  the  amplitude  of  the  secondary  lobes  for  the  full 
and  sparse  arrays  was  usually  less  than  10  dB. 

The  simulations  demonstrate  that  a  simple  linear  inter¬ 
polation  beamformer  can  produce  a  high  quality  radiation 
pattern  with  secondary  lobes  below  —60  dB  for  a  sparse 
array.  Such  a  system  is  feasible  with  present  technology 
and  could  theoretically  achieve  frame  rates  as  high  as  one 
thousand  2-D  images  per  second. 
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3.  DESIGN  AND  SIMULATION 


3.1  Circuit  Architecture 

The  first  step  in  designing  the  front  end  circuitry  for  the  system  is  determining 
the  most  difficult  and  most  important  specifications  to  achieve  and  select  parts  and 
topologies  which  will  make  these  specifications  easier  to  meet.  The  most  important 
specification  was  the  time-gain  control  with  45-dB  range.  The  second  most  important 
consideration  was  the  noise  requirement.  With  these  two  considerations  in  mind,  the 
Analog  Devices  part  AD604  was  selected  for  its  programmable  gain  from  0  dB  to  48 
dB  and  its  “ultralow”  input  referred  noise  (0.8  nV/VHz )  [8], 

The  third  most  important  design  consideration  was  the  high  voltage  protection 
circuitry.  A  topology  was  found  in  a  paper  by  Lockwood  and  Hazard  [9]  describing  a 
method  of  protecting  the  pre-amplifier  from  high-voltage  pulses.  The  final 
consideration  was  ensuring  that  the  analog-to-digital  converter’s  full  scale  range  was 
utilized.  A  simple  operational  amplifier  circuit  was  determined  to  have  an  adequate 
gain-bandwidth  for  this  application. 


High-Voltage  Pulses 


A/D 

Converter 


Transducer  Time  Gain  Control 


Figure  15.  Circuit  Architecture 
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3.2  Transducer  Simulation  Model 


In  order  to  optimally  design  the  high-voltage  protection  circuitry,  the  behavior 
of  the  ultrasonic  transducer  had  to  be  considered.  The  model  used  in  this  work, 
provided  by  Dr.  Lockwood,  is  based  on  the  KLM  model  for  ultrasonic  transducers 
[10,11],  The  model  uses  the  ABCD  network  parameters  of  basic  elements  to  build 
the  transducer.  These  ABCD  network  parameters  are  a  matrix  representation  of  a  two 
port  element.  These  two-by-two  matrices  relate  the  voltage  and  current  at  the  input  to 
the  voltage  and  current  at  the  output  by  the  equation 

Vm  _  A  B  Vom 

J IN  _  C  D  —  l OUT  _ 


An  entire  system  can  be  represented  by  multiplying  these  network  parameter  matrices 
together.  As  shown  in  Figure  16,  the  transducer  model  is  represented  by  a  series  bulk 
electrical  impedance  followed  by  a  transformer  from  electrical  to  acoustic  terms 
followed  by  an  acoustic  transmission  line  in  parallel  terminated  by  an  impedance 
representing  the  sound  transmitted  though  the  back  of  the  transducer  followed  by  two 
series  transmission  lines  representing  the  matching  layers  of  the  transducer.  The 
Matlab  code  used  for  simulation  can  be  found  in  Appendix  A. 

The  input  and  output  impedances,  as  well  as  the  maximum  output  voltage  of 
the  transducer,  are  critical  parameters.  The  input  and  output  impedances  of  the 
transducer  can  be  determined  from  the  product  of  all  of  the  ABCD  networks  by 
dividing  A  by  C  of  the  transmit  simulation  for  the  input  impedance  of  the  transducer, 
and  dividing  D  by  B  to  determine  the  output  impedance  [11].  The  input  and  output 
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impedance  are  the  same  if  the  transducer  medium  acts  like  an  infinitely  long 


transmission  line. 


Figure  16.  Transducer  KLM  Model 
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Figure  17.  Transducer  Input  and  Output  Impedance  vs.  Frequency 
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After  the  input  and  output  impedance  was  determined  in  Matlab,  a  SPICE 
model  was  constructed.  In  SPICE  the  transducer  was  modeled  as  a  transconductance 
with  magnitude  and  phase  from  a  lookup  table.  The  values  in  the  table  came  from  the 
results  of  the  Matlab  simulation.  The  input  and  output  impedance  versus  frequency  is 
plotted  in  Figure  17.  For  frequencies  less  than  3  MHz,  the  transducer  behaves  mostly 
like  a  capacitor  and  can  be  modeled  accurately  as  a  0.34-nF  capacitor  in  series  with  a 
40-Q  resistance. 

3.3  High-Voltage  Protection  Circuit 

The  topology  for  the  high-voltage  protection  circuit  described  by  Lockwood 
and  Hazard  [9]  uses  four  diodes  and  a  capacitor  to  isolate  the  high-voltage  pulse  from 
the  amplifiers  as  shown  in  Figure  18.  Pairs  of  diodes  are  connected  in  parallel  with 
opposite  polarities  to  form  two  bi-directional  voltage  clamps.  A  capacitor  couples  the 
output  from  the  transducer  to  the  amplifier  input.  One  diode  clamp  (D3/D4)  limits 
the  voltage  applied  to  the  pre-amplifier  during  the  high-voltage  pulsing,  while  the 
other  (D1/D2)  isolates  the  pulse  generator  from  the  relatively  weak  signals  returning 
from  the  transducer. 

The  circuit  operates  on  the  principle  that  a  high-voltage  pulse  turns  on  all  of 
the  diode  clamps,  so  there  is  a  conducting  path  from  the  pulse  generator  to  the 
transducer  and  a  path  from  the  amplifier  input  to  ground.  This  ensures  that  the 
voltage  at  the  input  of  the  amplifier  is  limited  to  the  on-voltage  of  the  diodes.  When 
the  high-voltage  pulse  is  released,  the  diode  clamps  turn  off  and  the  pulse  generator 
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input  does  not  have  a  conducting  path  to  the  circuit  and  the  transducer  is  capacitively 
coupled  to  the  input  of  the  amplifier. 

The  selection  of  the  capacitor  value  is  a  trade-off  between  the  power  coupled 
to  the  transducer  during  a  high-voltage  pulse  and  the  attenuation  of  the  received 
signal  from  the  transducer.  As  the  capacitor  size  increases,  the  power  received  by  the 

transducer  decreases  but  the  attenuation  of  the  signal  from  the  transducer  decreases. 

\ 

Another  consideration  when  selecting  the  capacitor  size  is  the  peak  current 
through  the  diodes.  An  order  of  magnitude  calculation  was  performed  to  determine  a 
relative  size  of  capacitor  required  for  the  circuit  to  work  well.  For  50-Q  coaxial 
cables  much  shorter  than  the  wavelength  of  the  signals  on  the  cable,  the  cables  can  be 
modeled  by  a  capacitor  equal  to  30  pF/ft  [14],  By  modeling  the  transducer  with  the 
SPICE  model  derived  from  the  Matlab  simulation,  modeling  the  coaxial  cables  as 
capacitors  of  value  300  pF,  and  using  the  input  impedance  of  the  AD604  (300  Kfl  in 
parallel  with  8  pF),  the  coupling  capacitor  needed  to  be  greater  than  500  pF.  Using 
500  pF  and  the  default  SPICE  diode  model,  it  was  determined  that  the  diodes  need  to 
carry  over  50  Amps  of  peak  current  during  a  high-voltage  pulse.  From  this 
consideration,  the  Vishay  Lite-on  S3  A  diodes  were  selected.  These  were  the  smallest 
100-A  maximum  surge  current  diodes  available. 

After  these  diodes  were  selected  simulations  were  run  to  determine  the 
optimum  capacitance  value  for  the  circuit,  so  that  the  current  through  the  D3/D4 
clamp  is  as  small  as  possible  during  a  high-voltage  pulse  and  the  voltage  at  the  input 
to  the  amplifier  has  the  least  attenuation  during  the  receive  phase.  The  value  selected 
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was  3  nF.  After  an  initial  prototype  was  constructed,  it  was  found  that  the  coupling 
capacitor  to  the  input  of  the  AD604  created  a  floating  node,  so  a  resistor  was  needed 
to  center  the  node  at  0  V.  The  resistor  size  selected  was  50  ft  This  value  was  chosen 
for  several  reasons.  The  first  reason  was  it’s  low  Johnson  thermal  noise,  given  by 

v2 

=  4km . 

A/ 

So,  a  50-ft  resistor  has  0.89  nV/  V Hz  of  thermal  noise.  The  second  reason  for 
choosing  50  ft  is  the  high  pass  filter  it  creates  with  the  3-nF  capacitor.  The  high  pass 
comer  frequency  for  this  is  approximately  1  MHz.  The  third  reason  for  selecting 
50  ft  was  to  terminate  the  50-ft  coaxial  cable  from  the  transducer  at  high  frequencies. 

After  the  addition  of  the  50-ft  resistor,  a  second  set  of  simulations  determined  that 
3  nF  was  still  the  optimum  value. 

Figure  20  shows  the  current  distribution  when  a  high-voltage  pulse  is  applied 
to  the  high-voltage  input.  Figures  19  and  21  show  the  main  function  of  the  circuit. 
Figure  19  shows  a  high-voltage  pulse  being  applied  to  the  transducer.  The  high- 
voltage  pulses  used  in  simulations  are  single-cycle  sine  waves  multiplied  by  raised- 
cosine  window  functions.  Figure  21  shows  the  input  of  the  amplifier  being  limited  to 
by  the  D3/D4  clamp  ±1  V  when  a  600-V  pulse  is  applied  to  the  transducer.  Figure  22 

shows  the  attenuation  of  a  signal  received  from  the  transducer  at  the  input  of  the 
amplifier. 
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Figure  18.  High-Voltage  Protection  Circuit  Diagram 


Figure  19.  Simulated  Voltage  Across  the  Transducer  When  High-Voltage  Pulse 

is  Applied 
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Figure  20 

.  Current  Through  D1/D2  Clamp  and  Cl  When  High-Voltage  Pulse  is 
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Figure  21.  Voltage  at  Input  of  Amplifier  When  a  High-Voltage  Pulse  is  Applied 
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Figure  22.  Voltage  at  Output  of  Transducer  Compared  to  Voltage  at  Input  of 
Amplifier  when  a  Received  Signal  is  Simulated 

3.4  AD604  Peripheral  Components 

The  Analog  Devices  part  AD604  [8]  is  a  24-pin  chip  with  two  ultra-low  noise, 
variable-gain  amplifiers.  Each  amplifier  has  1 1  pins  of  its  own  and  there  are  two  pins 
that  are  shared  by  the  two  amplifiers.  The  basic  architecture  of  each  amplifier 
includes  an  op-amp  in  a  non-inverting  configuration  with  a  gain  of  either  +14  dB  or 
+20  dB  depending  on  whether  the  feedback  pin  (FBK)  is  shorted  to  the  pre-amp 
output  pin  (PAO)  or  left  unconnected.  The  pre-amp  is  followed  by  a  differential 
attenuating  resistive  network  with  inputs  +DSX  and  -DSX.  This  network  is  internally 
connected  to  an  amplifier  with  a  gain  of  +34.4  dB  which  is  connected  to  the  output 
pin  (OUT).  The  attenuation  of  the  resistive  network  is  set  by  the  gain  control  pin 


(VGN).  Each  amplifier  has  three  power  supply  pins,  consisting  of  a  positive  supply 
(POS),  a  negative  supply  (NEG),  and  a  ground  (COM).  The  shared  pins  are  the 
common  mode  output  voltage  level  pin  (VOCM)  and  the  gain  control  reference 
voltage  (VREF). 


Figure  23.  Simplified  Block  Diagram  of  Single  Channel  of  AD604  [8] 

For  this  application  it  was  determined  that  the  +14  dB  gain  would  be  used  by 
shorting  the  FBK  pin  to  the  PAO  pin  because  the  lower  gain  allows  for  a  larger  signal 
to  be  amplified  in  case  the  amplitude  from  the  transducer  was  greater  than  expected. 
The  +DSX  and  -DSX  inputs  are  capacitively  coupled  to  the  PAO  and  ground 
respectively.  The  +DSX  and  -DSX  inputs  both  have  an  input  resistance  of  175  £2,  so 
the  coupling  capacitor  sizes  must  be  chosen  so  the  high-pass  comer  created  by  the 
capacitors  and  resistance  is  less  than  2  MHz.  20-nF  capacitors  were  chosen  to  set  the 
comer  at  50  kHz,  so  that  the  input  would  not  be  affected  by  the  roll-off  in  the  signal 
bandwidth  and  so  the  amplifiers  would  not  distort  the  signal.  The  VOCM  has  an 
internal  reference  circuit  that  sets  the  common  mode  of  the  output  to  2.5  V  when  a 
capacitor  is  connected  between  the  VOCM  and  ground,  so  a  0.1 -pF  capacitor  was 
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used.  The  VREF  pin  does  not  have  this  internal  reference,  so  a  2.5-V  external 
reference  is  generated  by  a  ZETEX  ZRA250A03  micro-power  voltage  reference,  as 
shown  in  Figure  24.  The  output  of  the  AD604  is  centered  about  2.5  V  set  by  the 
VOCM  voltage,  but  for  the  Pentek  model  9158  A/D  converter  the  input  is  assumed  to 
be  centered  about  0  V.  To  remedy  this,  the  output  of  the  AD604  is  capacitively 
coupled  to  the  second  stage  amplifier.  The  gain  control  (VGN)  is  provided  by  an 
external  power  supply  or  signal  generator. 

«VR 


*Gnd 

Figure  24.  ZRA250A03  Circuit  Diagram  and  Package  [12] 

3.5  Second  Stage  Amplifier  Design 

In  order  to  fill  the  full  scale  range  of  the  A/D  converter,  i.e.  -1  V  to  +1  V  a 
gain  of  at  least  +20  dB  was  required.  Since  the  polarity  of  the  signal  in  this 
application  does  not  matter,  and  inverting  amplifiers  have  superior  performance  to 
non-inverting,  an  inverting  op-amp  configuration  was  used,  as  shown  in  Figure  25. 
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Since  the  output  of  the  AD604  is  optimally  terminated  by  a  500-Q  resistor  the 
resistance  values  used  were  500  Q  and  5000  Q  The  coupling  capacitor  was  selected 
to  be  3  nF  between  the  output  of  the  AD604  and  the  500-0  resistor.  This  sets  the 
high-pass  comer  at  100  kHz.  The  gain-bandwidth  of  the  op-amp  is  required  to  be 
over  60  MHz,  so  at  a  gain  of  lOx,  the  bandwidth  of  the  amplifier  is  at  least  6  MHz. 
The  National  Semiconductor  CLC428  was  selected  for  its  low  noise  characteristics, 
160-MHz  gain-bandwidth,  and  the  fact  that  there  are  two  amplifiers  per  chip.  A  50-0 
resistor  was  added  to  the  output  of  the  amplifier  to  isolate  it  from  the  large 
capacitance  of  the  coaxial  cable  to  the  A/D  converter.  This  value  was  selected 
because  it  provides  proper  termination  for  the  cable. 


R2  5000 


Figure  25.  Inverting  Op-amp  Configuration. 
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3.6  Noise  Requirements 

The  maximum  input  to  the  A/D  converter  is  1  V.  The  maximum  noise  at  the 
output  to  the  amplifier  is  therefore  55  dB  less  than  1  V  or  1.7  mV.  The  maximum 
noise  at  the  input  of  the  amplifier  when  it  is  at  its  highest  gain  setting  is  65  dB  less 
than  2.5  mV,  or  1  pV.  Since  there  are  1 1  transducers  creating  each  transmit  pulse 

the  input  referred  noise  can  be  VTf  more  than  1  pV,  or  3.4  pV . 

3.7  Entire  Circuit  Simulations 

All  simulations  were  done  using  the  SPICE  models  of  the  parts  available  on 
the  manufacturers  web-sites.  Figure  26  shows  the  circuit  model  for  the  entire  circuit. 
The  magnitude  response  of  the  interface  circuit  from  the  transducer  input  to  the 
output  with  the  high-voltage  input  grounded  is  shown  in  Figure  27.  The  phase 
response  is  shown  in  Figure  28.  Figures  29  and  30  show  the  magnitude  and  phase 
response  for  the  signal  bandwidth.  By  taking  the  negative  derivative  of  the  phase  the 
group  delay  was  found  for  Figure  3 1 .  Figure  32  shows  the  input  referred  noise  of  the 
circuit  for  the  highest  gain  setting  verses  frequency.  The  total  input  referred  noise 
was  found  by  integrating  the  square  over  the  signal  bandwidth  and  taking  the  square 
root.  This  value  is  1.5  pVrms.  The  magnitude  response  and  phase  response  were 
zoomed  in  to  the  2  to  4  MHz  signal  band  to  show  the  flatness  of  the  magnitude 
response  and  linearity  of  the  phase  response.  Figure  33  shows  the  response  of  the 
circuit  when  a  high-voltage  is  applied  and  a  receive  signal  is  simulated  within  5  ps  of 
the  transmit  pulse. 
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Figure  28.  Phase  Response  of  Circuit 


Figure  29.  Magnitude  Response  in  2-4  MHz  Signal  Bandwidth 
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Figure  33.  Output  of  the  Circuit  when  a  High-Voltage  Pulse  is  Applied 
Followed  by  a  Receive  Signal 
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4.  FABRICATION  AND  TEST 


4.1  Prototype  Design  and  Fabrication 

A  two-channel  prototype  of  the  front-end  circuitry  was  designed  to  verify  the 
simulations.  Since  both  the  AD604  and  CLC428  are  dual  chips,  a  two-channel  design 
was  the  logical  choice,  but  requires  careful  consideration  of  cross-talk.  When 
designing  a  system  in  discrete  components  as  with  any  real  world  electrical  design, 
electromagnetic  interference  and  signal  coupling  are  major  concerns.  By  using 
surface  mount  components,  there  is  less  inductance  in  the  leads  of  devices  which 
means  there  is  less  of  an  antenna  to  pick  up  radio  frequency  interference.  Using  a 
ground  plane  with  surface  mount  components  helps  reduce  both  interference  and 
coupling.  The  ground  plane  does  not  allow  radio  waves  to  couple  into  the  circuit 
through  loops  created  by  power  supplies  going  through  chips.  The  ground  plane  also 
reduces  the  fringe  fields  of  components,  so  the  components  do  not  interfere  with  each 
other  and  create  oscillations. 

Printed  circuit  board  (PCB)  trace  placement  is  extremely  important  for 
reducing  interference.  Traces  carrying  small  signals  must  stay  far  away  from  traces 
carrying  large  interference  signals  to  avoid  capacitive  coupling.  Making  sure  that 
traces  do  not  form  loops  and/or  ground  plane  islands  reduce  the  inductance  that  picks 
up  external  electromagnetic  radiation. 

In  this  design,  special  attention  was  given  to  power  supply  lines.  Because  of 
the  large  gains  used,  any  power  supply  noise  or  interference  would  be  amplified 
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tremendously  and  appear  in  the  output  signal.  Using  bypass  capacitors  for  each  chip 
filters  out  some  of  the  high  frequency  noise  going  into  the  chips.  Also,  it  is  extremely 
important  to  not  have  the  power  supply  lines  to  close  to  the  output  of  the  amplifiers. 
This  keeps  the  output  from  coupling  to  the  power  supply,  creating  a  feedback  path 
which  could  lead  to  instability  and  oscillations. 

Ground  connections  are  also  very  important.  It  is  very  easy  to  inadvertently 
make  loops  in  ground  connections.  Special  care  must  be  taken  to  ensure  that  ground 
connections  are  made  to  the  correct  ground.  Signal  grounds  connected  to  instrument 
grounds,  etc.  can  make  ground  loops  via  the  building  steel.  Ground  loops  can  be  very 
inductive,  and  the  AD604  is  very  sensitive  to  inductance  on  the  ground  node. 

Choosing  the  proper  technology  for  each  of  the  components  in  the  design 
involved  some  trade-offs.  For  all  of  the  resistors  in  the  design,  surface-mount  thick- 
film  chip  resistors  were  used.  These  were  selected  because  of  the  availability  of  sizes 
needed,  small  size,  and  low  cost. 

The  capacitor  in  the  high  voltage  protection  circuit  is  a  ceramic  disc  capacitor 
rated  at  2  kV.  There  are  no  surface  mount  capacitors  available  that  can  withstand  the 
voltages  used  in  this  application.  The  rest  of  the  capacitors  in  the  design  are  surface- 
mount  multi-layer  ceramic  chip  capacitors,  which  have  adequate  high  frequency 
performance.  Thin  film  capacitors  were  used  for  a  preliminary  prototype,  but  when 
soldering  them  to  the  PCB  they  were  very  easy  to  melt  and  some  times  did  not  have 
the  correct  capacitance  after  soldering.  Ceramic  is  much  more  resistive  to  heat,  and 
therefore  much  easier  to  solder  to  the  PCB. 
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Figure  35.  Solder  Side  PCB  Layout 
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The  24-pin  Small  Outline  Integrated  Circuit  (SOIC)  package  of  the  AD604 
variable-gain  amplifier  was  chosen  because  it  is  surface  mount,  relatively  small,  and 
readily  available.  The  CLC428  op-amp  is  packaged  in  an  8-pin  SOIC.  The  power 
supplies,  grounds,  and  gain  control  voltages  are  all  connected  to  the  board  using 
banana  jacks.  The  transducer  input,  high-voltage  input,  and  output  are  connected  to 
the  board  by  50-0  PC  mount  BNC  connectors. 

The  printed  circuit  board  was  designed  using  the  MicroSim  PCBoards 
program,  a  part  of  the  MicroSim  DesignLab  suite.  The  board  has  two  layers  with  a 
minimum  trace  spacing  of  8  mils  (8/1000  inch).  Traces  in  this  design  were  not 
required  to  be  very  dense,  so  a  trace  width  of  30  mils  was  used  for  all  traces.  The 
boards  were  fabricated  using  a  T-Tech  Quick  Circuit  milling  machine.  The  copper 
traces  were  milled  out  of  the  double-sided  copper  clad  boards,  and  the  ground  plane 
was  left  where  the  boards  were  not  milled.  The  boards  were  populated  using  a  PACE 
soldering  station  with  tools  for  surface-mount  soldering. 

4.2  Test  setup 

The  two-channel  prototype  was  tested  in  a  lab  environment.  A  Hewlett 
Packard  (HP)  E3631  triple-output  power  supply  was  used  to  provide  ±5  V,  ground, 
and  the  gain  control  voltage,  Vgn.  The  +15-V  output  was  set  to  +5  V,  the  -15-V 
output  was  set  to  -5  V,  and  the  +6  V  output  was  varied  to  set  the  gain.  An  HP  33 14A 
function  generator/arbitrary  waveform  generator  was  used  to  generate  sine  waves 
used  in  the  tests.  This  model  was  used  because  of  its  ability  to  generate  signals  from 
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Figure  36.  Test  Setup 

0.1  mVp  to  10  Vp.  All  measurements  were  done  using  an  HP  54602  oscilloscope 
with  an  HP  54657  HP-EB  measurement  and  storage  module  connected.  All  three  of 
these  instruments  were  connected  to  a  Pentium  PC  with  an  HP  82341  HP-IB  interface 
card. 

Banana  plugs  were  used  to  connect  the  power  supplies  and  gain  control 
voltages  to  the  power  supply.  BNC  cables  were  used  to  connect  the  signal  generator 
and  oscilloscope  to  the  circuit.  Since  the  output  was  affected  by  capacitive  loading,  a 
relatively  short  1-ft.  BNC  cable  was  used  to  connect  the  output  of  the  circuit  to  the 
oscilloscope.  The  length  of  the  other  two  BNC  cables  that  connect  from  the  function 
generator  to  the  input  of  the  circuit  and  from  the  input  of  the  circuit  to  the 
oscilloscope  did  not  affect  the  performance  of  the  circuit. 

4.3  Frequency  Response  Tests 

A  program  was  written  in  Visual  Basic  in  Excel  to  control  stimulus  and 
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collect  data  for  measurement  of  the  frequency  response.  The  source  code  is  located 
in  Appendix  B.  The  HP-IB  interface  routines  were  called  to  control  the  oscilloscope 
and  function  generator  to  find  the  magnitude  and  phase  response  of  the  circuit.  The 
macro  sets  a  frequency  on  the  function  generator,  then  the  it  auto-scales  the 
oscilloscope  so  the  input  and  output  signals  are  measurable  by  the  scope  at  which 
point  it  queries  the  oscilloscope  to  find  the  amplitude  of  the  input  and  output  signals 
and  the  phase  difference  between  channel  1  and  channel  2.  The  macro  then  repeats 
this  procedure  for  all  of  the  frequencies  to  be  tested. 

Before  taking  measurements,  the  amplitude  of  the  input  sine  wave  from  the 
signal  generator  was  adjusted  so  that  the  output  was  1  Vp-p  when  a  3 -MHz  sine  wave 
was  applied  to  ensure  that  the  output  did  not  saturate  the  CLC428  op-amp.  The 
program  was  run  for  five  different  values  of  Vg„,  the  gain  control  voltage.  The  results 
of  the  magnitude  measurements  are  shown  in  Figure  37  and  38.  Figure  38  shows  that 
the  variation  in  gain  over  the  2-4  MHz  signal  bandwidth  is  less  than  0.7  dB.  Figure 
39  shows  the  3-dB  bandwidth  for  several  gain  settings,  calculated  from  the  magnitude 
response  of  the  circuit  from  Figure  37. 

The  phase  response  for  when  Vgn=0.5  V  is  shown  in  Figure  40  and  41.  Figures 
39  and  40  show  the  magnitude  and  phase  response  of  the  circuit  showing  more  detail 
around  the  signal  bandwidth.  In  Figure  41,  a  logarithmic  curve  fit  was  found  to 
accurately  match  the  phase  response  in  the  signal  bandwidth.  From  this  regression  the 
estimated  group  delay  was  found  by  computing  the  negative  derivative  of  the  phase 
with  respect  to  frequency.  This  result  is  shown  in  Figure  42. 
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Figure  37.  Measured  Magnitude  Response  for  Several  Gain  Settings 
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Figure  38.  Measured  Magnitude  Response  for  2-4  MHz  Signal  Bandwidth 
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4.4  Gain  Control,  Noise,  and  Cross-Talk  Tests 

The  dependence  of  gain  on  the  gain  control  voltage  was  measured.  The  setup 
was  exactly  the  same  as  for  the  frequency  response.  The  function  generator  was  set 
to  a  5-mVp-p,  3-MHz  sine  wave  and  connected  to  the  transducer  input  of  the  circuit. 
Measurements  were  taken  as  the  gain  control  voltage  was  increased  from  0.25  to  2.27 
V.  At  the  higher  gain  settings  (greater  than  55  dB),  the  CLC428  op-amps  began  to 
saturate,  and  the  function  generator  amplitude  was  decreased  to  0.5  mVp-p.  The 
results  are  shown  in  Figure  43. 


Figure  43.  Gain  versus  Gain  Control  Voltage 


A  test  was  performed  to  determine  the  total  input  referred  noise  of  the  system. 
In  this  test,  all  of  the  inputs  were  tied  to  ground  and  the  gain  was  set  to  maximum  (65 
dB).  The  output  was  then  analyzed  using  the  FFT  function  of  the  oscilloscope.  The 
input  referred  noise  in  the  circuit  bandwidth  was  calculated  by  dividing  the  FFT 
output  by  the  bandwidth  of  each  point  and  the  gain  in  the  signal  band.  The  total  noise 
was  found  by  integrating  the  mean  square  values  over  the  circuit  bandwidth.  The 
total  input  referred  noise  was  found  to  be  2.3  pVrms.  The  FFT  is  shown  in  Figure  44. 


Figure  44.  FFT  Showing  Input  Noise  Spectrum 

The  cross-talk  between  the  two  channels  of  the  prototype  was  measured.  In 
this  test,  a  sine  wave  was  applied  to  the  transducer  input  of  one  channel  of  the  circuit, 
with  a  magnitude  such  that  the  output  was  2  Vp-p.  The  output  from  the  second 


channel  was  measured  and  compared  to  the  output  of  the  first  channel.  The  cross-talk 
was  then  measured  for  several  gain  settings  and  the  results  are  shown  in  Figure  45. 
These  measurements  are  very  close  to  the  specifications  for  cross-talk  in  the  AD604. 


Figure  45.  Measured  Cross-Talk 


4.5  High-Voltage  Protection  Tests 

The  high-voltage  protection  circuitry  was  not  accurately  tested  because  the 
function  generators  necessary  to  apply  high  voltages  were  unavailable.  With  a  10-V 
pulse  applied  to  the  high  voltage  input,  the  voltage  at  the  transducer  output  was  about 
9  V  and  at  the  input  of  the  amplifier  was  about  1  V.  This  result  was  similar  to  what 
was  expected.  A  qualitative  trial  at  the  Cleveland  Clinic,  using  a  real  transducer  and 


the  actual  high-voltage  generator,  was  performed  and  the  high-voltage  protection 
circuitry  appeared  to  work  properly. 

4.6  Discussion 

The  results  of  the  frequency  response  tests  performed  on  the  circuit  were 
remarkably  close  to  the  simulation  results.  Figures  46,  47,  and  48  show  a  comparison 
of  the  two  results.  At  lower  frequencies,  the  shape  of  the  magnitude  response  of  both 
simulation  and  test  are  very  much  the  same,  but  at  the  higher  frequencies  the  test 
results  fall  off  more  rapidly  than  the  simulation.  This  roll-off  can  be  explained  by 
parasitic  capacitance  in  the  feedback  path  of  the  op-amp  and  in  the  coaxial  cables 
both  of  which  create  low-pass  filters.  The  prediction  of  gain  for  a  given  level  of  gain 
control  voltage  is  very  good  at  higher  gain  settings,  with  some  inaccuracy  at  lower 
settings. 


Figure  46.  Simulated  Magnitude  Response  Compared  to  Measured  Response 
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Figure  47.  Simulated  Phases  Response  Compared  to  Measured  Response 
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Figure  48.  Simulated  Group  Delay  Compared  to  Measured  Response 
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The  phase  response  is  also  affected  by  parasitic  capacitances  and  agreement 
between  simulation  and  measurement  is  best  at  lower  frequencies.  The  measured 
group  delay  is  greater  than  the  simulations  primarily  because  of  the  coaxial  cables. 
Of  key  importance  is  the  variation  of  the  group  delay  over  the  signal  bandwidth.  The 
simulated  variation  is  about  20  ns,  while  the  measured  variation  is  about  30  ns.  This 
quantity  in  the  test  is  less  than  10  ns  more  than  in  simulation.  The  input  referred 
noise  in  the  signal  bandwidth  was  twice  as  much  as  the  simulation  which  is  still 
within  specifications. 

The  bandwidth  of  the  circuit  is  a  bit  larger  than  what  is  necessary  for  this 
application.  The  bandwidth  extends  from  1  MHz  to  8  MHz  at  the  low  gain  setting 
and  from  1  MHz  to  6  MHz  at  the  high  gain  setting.  This  large  bandwidth  allows  for  a 
flatter  response  in  the  signal  band,  but  it  also  allows  more  total  noise  in  the  output 
signal. 

The  gain  control  voltage  controls  the  gain  in  a  linear-in-dB  fashion  as 
suggested  by  the  AD604  data  sheet  [7],  but  there  is  a  slight  non-linearity  at  both  the 
low  and  high  end  of  the  gain  settings.  When  the  signal  that  controls  the  gain  is 
created,  this  non-linearity  will  need  to  be  taken  into  account. 
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INTRODUCTION: 

The  purpose  of  this  report  is  to  describe  the  major  technical  risks  associated  with  the 
completion  of  the  objectives  of  the  subcontract  referenced  above  and  to  outline  the 
technical  approach  proposed  by  Tetrad  to  ameliorate  these  risks.  The  objective  of  the 
program  is  to  produce  a  working  64-element  phased  array  that  can  be  steered 
mechanically  in  the  plane  perpendicular  to  the  array  scan  plane.  Requirements  are  laid 
out  in  the  specifications  attached  as  Appendix  A. 


The  main  objectives  of  the  probe  are  as  follows: 


•  The  acoustical  design  must  be  capable  of  producing  good  images  over  a  90  degree 
sector  angle  without  excessive  artifact. 

•  The  probe  must  have  a  life  time  consistent  with  completion  of  a  thorough  evaluation 
of  the  imaging  system  and  probe. 

•  The  phased  array  must  be  capable  of  achieving  a  high  level  of  transmitted  energy  per 
element  to  achieve  the  high  signal  to  noise  ratio  required  for  effective  imaging  with  a 
synthetic  aperture  that  uses  a  relatively  small  transmitting  aperture  producing  a 
diverging  beam. 

Each  objective  carries  a  degree  of  technical  risk  that  must  be  addressed  in  the  program. 

They  are  addressed  separately  in  the  sections  that  follow. 

ACOUSTIC  DESIGN  ISSUES: 


The  mechanical  steering  of  the  array  requires  that  it  be  housed  inside  a  fluid-filled 
enclosure.  This  type  of  enclosure  is  commonly  used  for  mechanically  steered  single 
element  and  annular  array  transducers.  The  acoustic  problem  is  substantially  more 
difficult  for  a  phased  array,  however,  because  of  the  need  to  steer  the  beam  through  the 
enclosure. 
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Failure  to  achieve  a  good  acoustic  match  can  result  in  loss  of  signal  level  at  steered 
angles  and  in  reflections  inside  the  probe  that  cloud  the  near  field  of  the  image. 

Tetrad  has  begun  analysis  of  the  acoustic  requirements  for  acoustic  enclosures  of  varying 
shapes.  The  two  shapes  under  most  active  consideration  are  cylindrical  and  a  complex 
shape  approximating  a  toroidal  section.  See  Fig.  1. 

The  axis  of  the  cylindrically  shaped  enclosure  is  ideally  coincident  with  the  axis  about 
which  the  array  moves  mechanically.  This  causes  the  array  to  be  positioned  in  at  the 
same  distance  and  with  the  same  effective  window  shape  regardless  of  the  angle.  This 
enclosure  has  the  advantage  that  it  is  easy  to  construct  and  allows  the  array  to  be  very 
close  to  the  surface  of  the  enclosure.  It  has  the  disadvantage  that  beam  must  be  able  to 
pass  through  at  angles  ranging  from  normal  to  more  than  45  degrees  from  normal  without 
losing  significant  energy  or  creating  excessive  reflections.  This  requires  either  that  the 
enclosure  be  very  thin  or  that  it  be  well  matched  in  both  sound  velocity  and  acoustic 
impedance. 

The  toroidal  section  would  have  one  radius  with  a  center  of  curvature  on  the  axis  about 
which  the  array  moves  mechanically  and  another  radius  centered  at  the  center  of  the 
phased  array.  This  creates  a  shape  for  which  the  rays  deviate  from  normal  by  the 
smallest  amount  possible.  The  major  advantage  of  this  shape  is  that  the  acoustic  property 
constraints  are  less  severe  if  the  enclosure  is  placed  a  reasonable  distance  from  the  array. 
One  disadvantage  is  that  the  enclosure  must  be  molded  in  this  more  complicated  shape. 
More  importantly,  the  enclosure  is  further  from  the  array  so  any  reverberation  that  does 
occur  may  produce  artifacts  deeper  in  the  image. 

The  enclosure  has  two  interfaces  that  must  be  considered.  The  first  is  between  the  fluid 
and  the  enclosure  membrane  and  the  second  between  the  enclosure  membrane  and  the 
body.  Since  most  rigid  membrane  materials  have  sound  velocities  that  are  higher  than 
most  fluids,  it  is  the  fluid  to  membrane  interface  that  is  susceptible  to  critical  angle 
problems.  Fig.  2  is  a  plot  of  the  critical  angle  for  longitudinal  waves  as  a  function  of  the 
ratio  of  the  membrane  sound  speed  to  the  fluid  sound  speed.  It  is  encouraging  that  the 
ratio  of  sound  velocities  must  exceed  l  .4  before  critical  angle  drops  below  45  degrees.  A 
ratio  of  1.4  is  probably  achievable  in  common  fluids  and  membranes. 

It  is  also  important  to  consider  how  much  reflection  there  is  from  the  interface  and  how 
much  energy  will  be  transmitted  through  the  membrane.  To  estimate  this  we  have 
calculated  the  reflection  and  transmission  coefficients  for  the  inner  interface  ignoring 
mode  conversion  in  the  solid  material.  While  this  is  not  absolutely  correct  it  will  give  us 
the  right  order  of  magnitude  for  evaluation.  The  results  for  a  velocity  ratio  of  1 .3  and  an 
impedance  ratio  of  1.4  are  shown  in  Fig.  3.  Note  that  even  with  these  ratios,  the 
transmission  loss  is  less  than  0.3  dB  and  the  reflection  coefficient  is  less  than  -15  dB  at 
all  angles.  This  implies  that  the  major  loss  of  signal  at  an  angle,  other  than  the 
acceptance  angle  of  the  array  itself,  will  result  from  the  attenuation  in  the  membrane. 

The  path  through  the  membrane  for  a  45 -degree  angle  of  incidence  is  larger  than  the  path 
for  normal  incidence. 
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Selection  of  fluid  is  also  important.  The  velocity  and  impedance  must  match  the 
membrane  reasonably  well;  the  attenuation  must  be  such  that  the  reverberations  die  off 
quickly  but  low  enough  that  signal  to  noise  ratio  is  not  needlessly  reduced.  Some  fluids 
also  are  more  inclined  to  develop  small  bubbles  than  others.  This  relates  to  the  change  in 
air  solubility  over  the  temperature  range  of  intended  use. 

Based  on  this  analysis  it  will  be  Tetrad’s  approach  to  build  an  initial  probe  enclosure  with 
a  cylindrical  shape  containing  a  simple  phased  array.  We  will  do  either  real-time,  or 
synthetic  aperture  imaging  with  this  array  to  measure  the  signal  drop  off  and  to  examine 
the  reverberation  level.  If  this  approach  can  be  made  to  work  well  enough  with 
commercially  available  materials,  it  will  be  adopted  as  the  approach  of  choice.  If  there  is 
any  question  on  performance  level,  Tetrad  will  discuss  the  performance  the  necessity  for 
the  development  of  the  compound  curvature  enclosure  with  Goeff  Lockwood.  The 
development  of  the  compound  enclosure  could  become  substantially  more  expensive  than 
assumed  in  the  original  budgeting  and  may  need  to  be  assessed  relative  to  other  desirable 
features. 

PROBE  LIFE  TIME 

Most  of  the  components  of  the  mechanical  system  and  the  array  have  been  tested  in 
existing  systems  and  found  to  be  reliable.  The  motor,  for  example,  should  be  capable  of 
very  high  cycle  counts.  There  are  two  exceptions  to  this.  The  first  is  that  existing 
mechanical  probes  do  not  have  64  elements  that  must  be  brought  out  of  the  array  through 
a  moving  interface.  The  highest  element  counts  in  current  mechanical  systems  come 
from  annular  arrays  that  have  10  or  less  elements.  It  is  possible  with  this  low  number  of 
elements  to  use  coax  cables  for  the  connections.  With  the  phased  array  it  will  be 
necessary  to  go  to  a  flex  circuit  interconnection.  The  wire  connections  on  both  ends  of 
the  flex  can  be  well  potted  such  that  the  only  possible  mechanism  of  failure  would  be  in 
the  flex  leads  themselves. 

Flexes  are  often  used  in  high  cycle  count  equipment  such  as  in  laser  printers.  None  the 
less,  flex  life  is  a  parameter  that  should  be  tested  early  in  the  program  to  assure  that  there 
are  no  unexpected  failures.  To  do  this,  we  plan  to  have  a  traditional  phased  array  built  up 
and  placed  in  a  probe  for  life  testing  during  the  first  six  months  of  the  program.  This 
should  establish  that  the  life  expectancy  is  sufficient  or  guide  the  redesign  of  a  flex  for 
the  final  delivery. 

The  second  exception  is  that  the  phased  array  is  not  normally  exposed  to  probe  fluids. 

The  phased  arrays  will  be  soaked  in  probe  fluid  early  in  the  program  to  provide  assurance 
that  there  is  no  degradation  of  materials  or  acoustic  performance  in  the  fluid  environment. 

PROVIDING  HIGH  TRANSMITTED  ENERGY  PER  ELEMENT 

The  proposed  synthetic  aperture  approach  uses  a  reduced  number  of  array  elements  to 
provide  a  diverging  transmit  beam.  Analysis  done  by  Goeff  Lockwood  indicated  that  the 
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signal  to  noise  ratio  could  be  maintained  using  the  proposed  approach  if  the  transmit 
voltage  could  be  increased  by  a  factor  of  5. 

Since  phased  array  elements  may  be  pulsed  with  200  volts  or  more  the  voltage  required 
for  the  synthetic  aperture  approach  would  be  1000  volts  or  more  if  the  output  of  the 
piezoelectric  transducer  were  linear  with  respect  to  transmit  voltage.  However,  the 
piezoelectric  output  does  not  increase  linearly  with  applied  E  field,  so  it  is  possible  that 
even  more  voltage  might  be  required.  Furthermore,  the  maximum  E-field  is  high  enough 
that  it  could  cause  depoling  in  the  piezoelectric  material  if  it  were  applied  in  a  direction 
opposed  to  the  direction  of  polarization.  Consequently,  the  pulse  must  be  unipolar. 

These  constraints  produce  a  significant  design  challenge  for  the  transmitter. 

Tetrad’s  initial  proposal  was  to  provide  a  piezoelectric  array  with  5  ceramic  layers 
connected  electrically  in  parallel  and  acoustically  in  series.  This  would  allow  the 
transmitter  design  to  be  similar  to  that  of  a  standard  phased  array.  Transmitting  on  5 
layers  connected  electrically  in  parallel  with  1/5  the  voltage  that  might  be  applied  to  a 
single  layer  material  does  not  reduce  the  E  field.  This  means  that  a  unipolar  pulse  will 
still  be  required  and  the  reduction  in  acoustic  output  due  to  nonlinearity  will  still  be  seen 
with  the  multilayer  device. 

Fig.  4  shows  the  relationship  between  E  field  and  strain  for  a  material  similar  to  that 
which  would  be  used  for  the  array.  The  upper  branch  can  be  used  to  approximate  the 
strain  that  will  be  achieved  for  various  voltage  levels.  Given  the  stack  dimensions  at  3.5 
MHz,  the  E  field  for  200  V  and  a  single  layer  is  about  4.4  KV/cm.  This  produces  a  strain 
of  about  0.042%.  If  600  V  is  applied  to  a  single  layer  (or  200  V  to  a  triple  layer)  the  field 
is  13.2  KV/cm.  This  results  in  a  strain  of  0.102%  or  an  increase  of  2.4  times  in  strain  for 
an  increase  of  3  times  in  E  field.  If  1000  V  is  applied  to  a  single  layer  (or  200  V  to  five 
parallel  layers)  the  field  is  21.9  KV/cm  and  the  strain  is  0.138%.  This  is  an  increase  of  a 
factor  of  3.3  in  strain  for  an  increase  of  5  in  E  field.  It  also  appears  unlikely  that  any  E 
field  will  yield  an  increase  of  a  factor  of  5  in  strain. 

This  is  a  cause  of  some  concern  for  an  approach  requiring  five  times  the  acoustic  output 
whether  it  be  attempted  by  application  of  a  high  voltage  or  use  of  five  layers. 
Furthermore,  one  probable  multilayer  ceramic  vendor  objected  to  the  use  of  five  layers, 
saying  that  three  layers  was  probably  the  most  that  could  be  reliably  and  repeatably 
made. 

The  KLM  transducer  model  was  run  to  evaluate  the  possibility  of  using  three  instead  of 
five  layers  in  the  array.  The  results  were  rather  encouraging.  The  model  predicts  that  the 
signal  to  noise  ratio  in  a  three  layer  array  is  16  dB  better  than  that  for  a  single  layer  array. 
When  the  strain  nonlinearity  is  taken  into  account  the  improvement  is  14  dB,  which  is 
exactly  the  factor  of  5  required  by  the  analysis  done  by  Goeff  Lockwood.  Much  of  this 
improvement  takes  place  on  receive  since  the  layered  elements  are  much  better  at  driving 
the  cable  than  are  the  small  single  layer  elements.  Figs.  5  and  6  are  the  input/results 
pages  from  the  KLM  model  for  the  single  layer  and  triple  layer  array  elements 
respectively. 
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There  is  still  some  risk  in  using  even  the  3 -layer  system  since  3 -layered  ceramic  is  still 
not  commonly  available  in  the  United  States.  The  primary  risk  is  that  the  material  will 
not  tolerate  dicing  without  damaging  at  least  a  few  elements.  To  reduce  that  risk  Tetrad 
plans  to  purchase  stacked  ceramic  from  multiple  sources  and  evaluate  the  dicability  of  the 
materials  within  the  first  eight  months  of  the  program.  If  good  results  can  be  established, 
work  with  the  multilayers  will  proceed.  If  the  results  are  not  good,  then  Tetrad  and 
Cleveland  Clinic  will  reassess  the  requirements  and  attempt  to  determine  a  combination 
of  electronics  and  arrays  that  will  meet  the  signal  to  noise  requirements.  The 
requirements  may  again  require  careful  financial  analysis  to  avoid  development  that  goes 
beyond  the  funding  available. 

At  the  time  of  this  report  Tetrad  has  received  the  first  set  of  three  layer  ceramic,  has  diced 
1 8  elements  without  apparent  damage.  The  performance  appears  to  be  good. 


Fig.  1:  Cylindrical  window  at  left.  Mechanical  pivot  is  on  the  cylindrical  axis.  Toroidal  section 
membrane  at  right.  Mechanical  pivot  is  at  the  center  of  the  toroid.  Phased  array  pivot  (array 
center)  moves  on  the  circle  that  describes  the  center  of  the  toroidal  tube. 
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The  main  efforts  of  the  program  during  this  quarter  had  to  do  with  subcontracting  the 
mechanical  mechanism  and  continuing  to  investigate  multilayer  phased  arrays. 

Regarding  the  mechanism;  the  specifications  were  completed  and  the  subcontract  signed 
in  early  November.  A  single  layer  phased  array  with  suitable  performance  for  testing  the 
probe  acoustics  has  been  developed  for  testing  the  acoustic  performance  of  the 
mechanical  scanner  dome  and  fluid  system.  The  flex  circuit  for  attaching  the  array  has 
been  designed  fabricated.  Acoustic  testing  of  membrane  materials  and  fluids  has  been 
completed  for  materials  that  should  be  compatible  the  use  of  a  phased  array  with  a 
cylindrical  dome.  Initial  prototypes  for  in-situ  acoustic  testing  are  being  constructed. 

A  mechanism  and  motor  driver  has  been  configured  for  achieving  20  frames  per  second 
and  life  testing  is  in  progress. 

Work  with  the  multilayer  phased  array  is  slightly  behind  the  original  schedule.  The 
single  layer  phased  array  that  was  constructed  as  a  baseline  showed  approximately  7  dB 
less  sensitivity  than  comparable  commercial  units.  The  cause  of  this  sensitivity  loss  has 
been  identified  and  single  units  with  sensitivity  comparable  to  other  commercial  units  can 
now  be  made  consistently.  Alternate  lensing  studies  are  underway  to  improve  the 
sensitivity  more. 

One  unit  has  been  constructed  using  multilayer  ceramic.  The  geometry  for  the  multilayer 
is  not  exactly  the  same  as  it  was  for  the  single  layer  units  but  the  increase  in  sensitivity  on 
the  best  elements  is  on  the  order  of  12  dB  more  sensitive  on  the  multilayer  units. 
Unfortunately  some  elements  showed  signs  of  losing  internal  electrodes.  A  second  unit  is 
being  built  using  more  care  in  selecting,  handling,  and  dicing  the  multilayer  units.  Some 
units  also  showed  signs  of  depoling.  This  appears  to  have  been  corrected  after  they  were 
returned  to  the  vendor  to  undergo  a  revised  poling  process. 

It  is  our  goal  to  reach  a  decision  on  the  use  of  multilayers  by  the  end  of  January,  2000. 
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Abstract-  A  synthetic  aperture  (SA) 
beamformer  is  theoretically  capable  of  collecting 
volume  data  sets  in  real-time.  Since  SA  imaging 
requires  phase  coherence  over  multiple  transmit 
events,  relative  motion  between  the  transducer  and 
target  is  a  concern,  especially  for  a  cardiac  imaging 
system.  The  sensitivity  to  motion  of  a  sparse  transmit 
SA  beamformer  was  evaluated  by  simulating  a  sector 
scan  of  a  moving  point  target.  The  results  were 
quantified  by  plotting  the  maximum  value  of  the 
beamformed  signal  as  function  of  scan  angle. 
Simulations  were  used  to  examine  motion  in  both  the 
axial  and  lateral  directions.  Lateral  motion,  for  the 
speeds  simulated  (<2  m/s),  did  not  significantly  change 
the  secondary  lobes  in  the  radiation  patterns.  Motion 
in  the  axial  direction  produced  more  distortion  of  the 
radiation  pattern;  however,  the  secondary  lobes  for 
targets  moving  at  speeds  representative  of  cardiac 
motion  (<80  mm/s)  were  not  significantly  different 
than  those  for  stationary  targets.  For  larger  speeds 
(400-800  mm/s)  the  secondary  lobes  rose  about  4  dB. 
These  findings  were  verified  experimentally  using  a 
2.6  MHz  array.  The  results  indicate  that  SA 
beamforming  with  a  sparse  transmit  array  may  be 
feasible  for  cardiac  imaging. 

I.  Introduction 

Real-time  3-D  ultrasound  imaging  provides 
unique  clinical  information.  We  have  previously 
described  a  method  for  real-time  3-D  imaging  using  a 
sparse  linear  phased  array  and  synthetic  aperture  (S A) 
beamforming  [1],  [2].  A  real-time  3D  data  set  is 
collected  by  mechanically  rocking  the  array  in  the 
elevation  direction  while  rapidly  acquiring  2D  images. 
Figure  1  depicts  the  collection  sequence  for  a  SA 
image.  A  single  diverging  pulse  is  transmitted  and  the 


Transmit  Receive 

Figure  1:  Collection  sequence  for  synthetic 

aperture  system 

ehcoes  are  simultaneously  collected  on  all  receive 
channels.  The  signal  for  each  channel  is  digitized  and 
stored.  This  procedure  is  repeated  for  each  transmit 
events.  The  stored  signals  are  then  delayed,  apodized, 
and  summed  digitally  to  form  all  of  the  lines  in  the 
image  simultaneously.  Provided  the  transducer  and 
target  tissue  are  stationary  during  data  collection,  the 
phase  coherence  necessary  for  proper  beamforming  is 
maintainted.  However,  if  the  target  or  transducer 
moves  during  image  acquisition,  the  radiation  pattern 
will  be  degraded.  Cardiac  imaging  is  one  application 
in  which  target  motion  is  present.  The  velocity  of  the 
heart  wall  has  been  reported  to  be  between  50  mm/s 
and  80  mm/s  in  systole  and  lower  in  diastole  [3].  Since 
the  transducer  is  continuously  rocked  during  the 
collection  of  a  3D  data  set,  transducer  motion  is  also  a 
concern. 

The  problem  of  motion  in  S A  imaging  systems 
has  been  investigated  previously,  and  several 
researchers  have  proposed  algorithms  to  restore  phase 
coherence  in  the  presence  of  tissue  motion  [4], [5].  In 
this  paper,  we  investigate  how  much  motion  can  be 
tolerated  by  a  sparse  transmit  SA  imaging  system 


before  the  loss  of  phase  coherence  significantly 
deteriorates  the  radiation  pattern.  The  motion 
sensitivity  was  investigated  by  simulating  a  sector 
scan  of  a  moving  point  target.  This  computer 
simulation  was  validated  using  an  experimental  rf  data 
set. 

II.  Experimental  Data  Collection 

A  96  element  transducer  array  (General 
Electric,  model  227S),  with  a  center  frequency  of  2.6 
MHz  and  0.28  mm  element  spacing,  was  used  to 
collect  the  experimental  rf  data  sets.  64  elements  of 
the  array  were  used.  The  array  was  positioned  at  the 
top  of  a  water  tank  using  a  precision  translation  stage 
which  allowed  for  precise  up  and  down  movement. 
The  target  was  a  small  sealed  plastic  cylinder, 
positioned  7  cm  (~f/4.3)  from  the  face  of  the  array 
with  the  long  axis  perpendicular  to  the  image  plane. 

Data  collection  followed  the  sequence  shown 
in  figure  1 .  A  single  element  was  excited  (+/-  40  V) 
and  the  echoes  were  recorded  on  each  of  the  64 
elements  individually.  (In  a  real-time  system  the 
signals  would  be  recorded  simultaneously.)  This 
procedure  was  repeated  for  each  of  three  transmit 
events.  The  signals  received  by  the  array  elements 
were  amplified  by  a  low  noise  pre-amplifier  (Miteq) 
and  digitized  at  50  MHz  and  8  bits  by  a  digital 
oscilloscope  (H.P.  54522A).  An  averaging  feature  of 
oscilloscope  was  used  to  increase  the  number  of 
effective  bits  to  10.  The  digital  data  was  transmitted 
to  a  PC,  where  the  dc  component  was  removed  before 
beamforming. 

To  simulate  the  effect  of  motion,  the  target 
was  moved  a  specified  distance  between  transmit 
events.  This  distance  was  calculated  by  multiplying 
the  velocity  to  be  simulated  by  the  time  between 
transmits  (200  /zs).  Experimental  rf  data  sets  were 
collected  for  a  stationary  target  and  targets  with  axial 
speeds  of  400  mm/s  and  800  mm/s. 

Collecting  an  experimental  rf  data  set  for 
synthetic  aperture  is  very  tedious.  To  facilitate 
exploring  of  a  wider  range  of  parameters,  a  computer 
was  used  to  generate  simulated  rf  data  sets.  While  this 
simulation  does  not  completely  model  the 
experimental  environment,  it  provides  useful 
information  for  predicting  the  performance  of  the 
system.  The  experimental  data  was  used  to  validate 
the  results  of  the  computer  simulations. 
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Figure  2:  Simulated  and  Experimental  pulse,  a) 
pulse  b)  spectrum. 

III.  Computer  Simulations  and  Beamforming 

We  simulated  an  array  with  64  elements 
spaced  at  0.28  mm,  and  the  center  frequency  of  the 
simulated  rf  pulses  was  2.6  MHz.  The  shape  of  the 
pulse  used  in  the  computer  simulations  was  matched 
to  the  experimental  pulse  so  that  the  resulting 
radiation  patterns  could  be  compared.  Figure  2 
compares  the  simulated  pulse  with  an  experimental 
pulse,  both  in  the  time  and  frequency  domains.  The 
same  sparse  transmit  aperture  used  in  the  experiment 
was  simulated.  The  signals  for  a  point  target  placed  at 
f/4.3  and  -2.2°  were  calculated  for  each  of  the  active 
transmitter-receiver  pairs  and  sampled  at  50  MHz  and 
10  bits.  These  signals  were  stored  and  used  as  an 
input  to  the  same  beamforming  program  that  was  used 
for  the  experimental  data. 

The  beamforming  program  generates  a  sector 
scan  from  -45°  to  45°.  Each  point  in  the  sector  scan 
was  generated  by  summing  appropriately  delayed 
signals  from  each  channel  for  each  transmit.  The 
delayed  signal  values  were  calculated  by  linearly 
interpolating  between  the  two  surrounding  digitized 
data  points  [1].  A  raised  cosine  apodization  was  used 
on  receive  and  the  outer  transmits  were  given  one  half 
the  weight  of  the  central  transmit.  Second  order 
sampling  was  used  to  approximate  the  envelope  of  the 
beamformed  signals.  A  radiation  pattern  was 
generated  by  plotting  the  maximum  value  of  the 
demodulated  signal  along  each  scan  line  as  a  function 
of  the  scan  angle. 


IV.  Results 

A.  Lateral  Motion 

A  rocking  array,  for  real-time  3D  imaging  at 
20  volumes  per  second,  will  produce  lateral  velocities 
of  about  1.6  m/s  at  the  depth  of  15  cm.  Figure  3 
shows  the  effect  lateral  motion  on  the  main  lobe  of  the 
radiation  pattern  for  velocities  of  0,  1  m/s,  and  2  m/s. 
The  main  lobe  is  shifted  and  slightly  wider  in  the 
presence  of  lateral  motion.  The  secondary  lobes 
remain  unchanged  in  the  presence  of  the  lateral 
motion. 


Figure  3:  Effects  of  simulated  lateral  motion  on  the 
radiation  pattern. 
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Figure  4:  Comparison  of  the  radiation  patterns  for 
axial  velocity  of  0  mm/s  and  100  mm/s. 


B.  Axial  Motion 

Axial  motion  causes  a  more  serious  disruption 
of  the  phase  coherence  required  for  beamforming. 
This  can  result  in  a  deterioration  of  the  resolution  and 


contrast  in  a  SA  image.  Simulations  of  the  radiation 
pattern  for  a  point  target  moving  with  an  axial  velocity 
of  0  mm/s  and  100  mm/s  (representative  of  cardiac 
motion)  are  shown  in  figure  4.  The  results  show  no 
significant  difference  in  the  level  of  the  secondary 
lobes  only  a  slight  distortion  of  the  main  beam. 

Figure  5  compares  the  theoretical  and 
experimental  radiation  patterns  for  stationary  targets 
and  target  velocities  of  400  mm/s  and  800mm/s. 
There  is  good  agreement  between  the  theory  and 
experiment  down  to  about  -50  dB.  The  secondary 
lobes  are  below  -54  dB  for  all  cases.  The  unique 
structures  in  the  main  lobe,  caused  by  the  axial 
motion,  are  seen  in  both  the  theory  and  experiment, 
though  not  quite  as  pronounced  in  the  experiment. 
The  secondary  lobes  for  the  theory  are  about  6  dB 
lower  than  the  experimental  data.  Since  the 
experimental  data  was  collected  over  a  long  period  of 
time,  mechanical  vibrations  and  thermal  fluctuations 
may  have  contributed  to  this  difference.  Axial 

motion  causes  a  shift  in  the  location  of  the  peak  and 
slight  widening  of  the  main  beam,  in  addition  to  the 
changes  in  the  fine  structure  of  the  main  beam.  The 
level  of  the  secondary  lobes  for  a  stationary  target  is 
about  -64  dB  for  the  simulation  and  -58  dB  for  the 
experiment.  The  secondary  lobes  for  cases  with  axial 
motion  of  400  mm/s  and  800  mm/s  are  about  4  dB 
higher  than  those  for  the  stationary  target.  However, 
the  level  of  the  secondary  lobes  for  an  axial  target 
velocity  of  400  mm/s  is  not  significantly  different 
from  level  for  an  axial  velocity  of  800  mm/s.  This  is 
true  for  both  theory  and  experiment.  The  level  of  the 
secondary  lobes  in  the  presence  of  axial  motion  is  -60 
dB  for  theory  and  -54  dB  for  experiment. 

V.  Summary  and  Conclusions 

In  this  paper,  we  have  presented  experimental 
and  simulated  evidence  that  motion  can  be  tolerated  in 
a  sparse  transmit  synthetic  aperture  imaging  system. 
In  both  experiment  and  simulation  a  2.6  MHz,  64 
element  array  with  signals  digitized  at  10  bits  and  50 
MHz  was  evaluated  with  a  point  source  at  f/4.3. 
Linear  interpolation  was  used  to  increase  the  delay 
accuracy  in  beamforming. 

The  results  of  the  simulations  show  that 
lateral  motion  does  not  significantly  distort  the 
radiation  pattern  for  the  velocities  simulated  (<  2m/s). 
There  is  no  rise  in  the  level  of  the  secondary  lobes 
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Figure  5:  Effects  of  axial  motion  on  the  radiation 
pattern.  Theory  and  experiment  for  a)  no  motion 
b)  400  mm/s  and  c)  800  mm/s. 


and  only  a  minor  shift  and  widening  of  the  main  beam. 

While  axial  motion  is  a  more  serious  concern, 
simulations  reveal  that  at  speeds  typical  of  cardiac 
motion  (80  mm/s)  there  is  little  distortion  of  the  main 
beam  and  virtually  no  increase  in  the  level  of  the 
secondary  lobes.  At  higher  speeds  (400  mms-800 
mm/s),  there  is  about  a  4  dB  increase  in  the  level  of 
the  secondary  lobes.  Axial  motion  also  produces  fine 
structure  in  the  main  lobe  and  widens  the  beam 
slightly.  These  results  were  seen  in  both  experimental 
data  and  simulation.  The  results  suggests  that  real¬ 
time  3D  imaging  using  a  synthetic  aperture 
beamformer  is  feasible  for  cardiac  imaging. 
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ABSTRACT 

Real-time  volume  ultrasound  imaging  of  the  heart  is  a  new  trend,  and  the  registration  of  acquired  volume  framesets  is  clinically 
important.  This  registration  may  be  accomplished  by  processing  a  selected  pair  of  volume  frames  having  identical  cardiac  phase 
(preferably  end  diastolic)  from  two  framesets.  The  registration  solves  for  the  optimal  rigid  transformation  between  selected 
volumes  through  maximization  of  mutual  information,  a  voxel  similarity  measure.  The  accuracy  of  registration  was  estimated 
through  retrieval  of  an  artificially  introduced  misalignment.  Two  volume  frames,  belonging  to  the  same  frameset  and  separated 
in  time  by  250  ms,  were  selected.  The  secondary  volume  was  translated  by  seven  voxels  along  each  axis  and  rotated  by  seven 
degrees  about  each  axis  relative  to  the  primary  prior  to  registration.  The  translational  mismatch  upon  registration  was  within  one 
voxel  and  the  rotational  mismatch  less  than  two  degrees.  Reduction  of  the  speckle  noise  by  spatio-temporal  averaging  followed 
by  intensity  binning  was  a  key  step  in  successful  application  of  mutual  information  approach  to  ultrasound  imaging.  The 
application  of  our  method  to  nine  framesets  arising  from  four  different  patients  demonstrates  the  feasibility  of  using  of  mutual 
information  for  automatic  registration  of  cardiac  ultrasound  data. 

Keywords:  3D  Registration,  3D  Ultrasound,  Mutual  Information. 

1.  INTRODUCTION 

Real  time  three  dimensional  (3D)  acquisition  is  one  of  the  latest  developments  in  the  field  of  ultrasound  imaging.  The  traditional 
ultrasound  transducers  produce  planar  images,  typically  at  20-30  frames  per  second.  Before  the  advent  of  real  time  3D  ultrasound 
imaging,  attempts  were  made  to  create  3D  images  through  stacking  of  planar  images  by  carefully  sweeping  a  traditional 
ultrasound  transducer  over  a  3D  region1.  This  approach  was  however  prone  to  patient  movement  and  transducer  localization 
errors,  and  more  importantly  is  not  fast  enough  to  image  a  moving  anatomy  such  as  the  heart  in  three  dimensions.  In  contrast, 
real  time  3D  ultrasound  acquisition  is  capable  of  producing  approximately  20-30  volumes  per  second  and  truly  capturing  the 
motion  of  the  heart  with  a  very  high  temporal  resolution.  In  essence,  real-time  3D  ultrasound  imaging  is  four  dimensional  (4D) 
imaging  with  time  being  the  fourth  dimension. 

Real  time  3D  ultrasound  enables  many  new  applications;  one  of  which  is  registration  of  4D  images  or  volume  framesets. 
Registration  of  two  ultrasound  volume  framesets  taken  before  and  after  a  clinically  significant  event  can  provide  important 
diagnostic  information.  In  this  paper  we  present  our  ongoing  work  on  the  registration  of  cardiac  volume  framesets.  The  eventual 
goal  of  the  present  work  is  to  visualize  the  deviations  in  cardiac  size  and  shape  between  two  ultrasound  examinations.  The 
registration  is  therefore  limited  to  rigid  geometric  transformation  described  by  six  parameters  -  three  associated  with  translation 
and  three  with  rotation.  An  additional  registration  parameter  is  the  volume  frame  number  that  depends  on  the  accompanying 
electrocardiogram  (ECG)  and  is  described  in  detail  in  the  next  section. 

There  are  two  main  approaches  to  image  registration:  segmentation  based  and  voxel  similarity  based2'3.  Segmentation  based 
image  registration  seeks  to  first  identify  matching  structures  in  the  two  datasets  and  then  fit  the  identified  matching  structures 
using  a  fitting  criterion.  Segmentation  of  real  time  3D  ultrasound  images  is  a  complex  problem  in  its  own  right  and  the  degree 
of  true  correspondence  between  the  matching  structures  is  always  of  suspect.  This  approach  was  therefore  considered  unsuitable 
for  the  present  registration  application.  The  second  approach  to  image  registration  exploits  the  underlying  data  (pixel  or  voxel) 
similarity  and  does  not  depend  on  prior  segmentation.  This  approach  is  also  amenable  to  both  intra-  and  inter-modality  image 
registration  with  full  automation.  A  data  similarity  approach  was  deemed  most  fitting  for  the  problem  at  hand. 
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Among  many  data  similarity  measures,  mutual  information  (MI)  has  been  demonstrated  to  be  most  effective  and  robust  for  multi¬ 
modality  registration  of  a  structural  imaging  modality  (MRI  and  CT)  with  a  functional  imaging  modality  (PET,  SPECT,  etc.)4*5. 
Our  choice  of  mutual  information  as  the  measure  of  data  similarity  was  based  on  its  reported  success.  We  do  note  that  we  have 
found  only  one  example  of  the  use  of  mutual  information  in  ultrasound  image  registration6  so  far.  In  general,  image  registration 
involving  ultrasound  data  has  not  been  investigated  actively.  The  present  study  is  therefore  an  investigation  into  a)  the  feasibility 
of  using  mutual  information  for  ultrasound  image  registration  and  b)  identification  of  a  procedure  to  compute  mutual  information 
for  most  effective  results. 


2.  REGISTRATION  APPROACH 

Before  describing  our  technique,  a  word  on  terminology  is  essential  due  to  the  new  imaging  modality  involved.  We  will  refer 
real  time  3D  ultrasound  imaging  simply  as  3D  ultrasound  imaging  in  the  remainder  of  the  paper.  3D  ultrasound  image  acquisition 
produces  a  sequence  of  volume  frames  which  will  be  called  a  frameset.  A  frameset  will  be  considered  synonymous  with  a  4D 
image.  A  volume  frame  within  a  frameset  will  be  referred  to  as  a  frame. 

Our  objective  has  been  to  register  a  pair  of  cardiac  ultrasound  framesets  or  in  other  words  perform  4D  registration.  We  assume 
that  a  single  3D  rigid  transformation  differentiates  the  two  framesets.  This  is  a  safe  and  reasonable  assumption  given  the  short 
image  acquisition  duration  to  allow  any  fluctuation  in  the  transducer  orientation  and  breathing-related  movement  artifacts.  Owing 
to  rapid  acquisition  rate  (20-30  frames  per  second)  of  3D  ultrasound  imaging,  a  typical  examination  lasts  less  than  two  seconds 
giving  forth  more  than  one  full  cardiac  cycle  of  data.  The  4D  registration  problem  can  therefore  be  seen  as  consisting  of  two 
distinct  parts:  temporal  registration  and  3D  spatial  registration. 

Temporal  registration  -  Temporal  registration  is  the  process  of  establishing  frame-by-frame  correspondence  based  on  cardiac 
phase  between  two  given  framesets.  This  step  would  not  have  been  required  if  the  acquisitions  were  ECG  triggered.  In  the 
present  study  however,  temporal  registration  ensured  that  a  frame  showing  for  example  the  onset  of  diastole  in  the  first 
frameset  corresponded  with  a  frame  with  the  same  phase  in  the  second  sequence.  When  the  frame  rate  is  identical  between 
the  two  framesets,  a  one-point  synchronization  suffices  frame-by-frame  matching  for  all  the  remaining  frames.  If  the  frame 
rate  is  different,  data  interpolation  may  be  required  along  the  time  axis.  The  frame  rate  was  kept  constant  in  the  present  study 
eliminating  any  need  for  data  interpolation  for  temporal  registration. 

3D  spatial  registration  -  3D  spatial  registration  is  the  search  for  the  optimal  3D  rigid  transformation  between  a  pair  of 
"temporally  registered"  frames,  called  the  primary  and  the  secondary.  This  rigid  transformation  when  applied  to  the 
secondary  puts  it  in  the  reference  frame  of  the  primary.  As  explained  above,  we  assume  that  this  same  3D  rigid 
transformation  describes  the  misalignment  between  each  pair  of  frames. 

Any  cardiac  phase  can  be  chosen  to  perform  temporal  registration  and  the  subsequent  3D  spatial  registration.  Theoretically,  the 
frame  rate  of  3D  ultrasound  imaging  is  high  enough  to  capture  short-lived  cardiac  events  such  the  onset  of  systole,  and  all  phases 
should  yield  the  same  result.  Practically,  the  choice  of  end  diastolic  phase  wherein  the  heart  is  typically  stationary  in  human  for 
200-250  ms  is  the  most  prudent  selection.  A  small  time  shift  error  in  temporal  registration  will  have  little  influence  on  the 
accuracy  of  3D  spatial  registration.  Furthermore,  since  the  temporal  resolution  of  3D  ultrasound  imaging  is  not  infinite,  the  end 
diastolic  images  are  expected  to  have  the  least  motion  artifact.  The  desired  central  end  diastolic  frame  can  be  determined  by 
analyzing  ECG  that  accompanies  routine  cardiac  echocardiographic  examinations.  We  also  present  an  alternate  approach  to 
determining  the  most  stationary  end  diastolic  phase  by  computing  mutual  information  between  adjacent  frames.  The  more 
stationary  the  heart,  the  greater  is  the  correlation,  and  hence,  the  higher  is  the  mutual  information  between  two  adjacent  frames. 
The  frame  with  the  highest  "adjacent  frame"  mutual  information  along  the  sequence  corresponds  well  with  the  central  end 
diastolic  phase  given  by  ECG. 

The  mutual  information  I(A,B)  between  two  volume  frames  A  and  B  depends  on  individual  probability  density  functions  p(a) 
and  p(b),  and  the  joint  probability  density  function  p(a,b)  of  voxel  intensities  in  the  overlapping  zone  of  A  and  B.  Mathematically, 
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Both  a  and  b  vary  in  the  interval  [0,M]  where  M  is  highest  possible  voxel  intensity  value.  The  intensity  distributions  are 
accumulated  as  below: 
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where  /*(«)  and  n(b)  are  number  of  A  and  B  voxels  with  intensity  a  and  b  respectively  in  the  overlapping  zone,  and  n(a,b)  is  a 
similar  quantity  for  the  joint  distribution.  N  is  the  total  number  of  voxels  in  the  overlapping  zone. 


The  general  approach  of  3D  spatial  registration  is  the  following.  The  primary  frame  of  the  selected  pair  is  kept  fixed  spatially 
whereas  the  secondary  is  subjected  to  a  series  of  rigid  transformations.  Each  successive  transformation  brings  the  secondary 
progressively  in  alignment  with  the  primary  according  to  a  criterion  function,  i.e.,  mutual  information  in  the  present  study.  The 
process  is  said  to  have  converged  when  mutual  information  is  maximized  and  any  further  transformation  yields  a  smaller  mutual 
information  value.  A  rigid  transformation  in  3D  is  described  by  a  matrix  of  six  independent  parameters  -  three  controlling 
translation  and  three  rotation.  During  successive  iterations,  original  voxels  of  the  secondary  do  not  generally  coincide  with  those 
of  the  primary.  The  calculation  of  the  criterion  function  is  preceded  by  interpolation  (trilinear  interpolation  in  the  present  study) 
of  the  transformed  secondary. 


The  3D  spatial  registration  involves  an  iterative  search  for  the  optimal  rigid  transformation  that  maximizes  MI  in  a  six¬ 
dimensional  parameter  space.  There  are  a  few  important  issues  to  keep  in  mind  when  selecting  the  search  algorithm.  It  has  been 
shown  that  the  hypersurface  defined  by  the  MI  function  is  generally  uneven  with  a  lot  of  ripples5.  Our  findings  confirm  this.  The 
search  algorithm  should  therefore  be  robust  to  local  maxima  and  converge  only  to  the  global  maximum.  A  robust  search 
algorithm  however  takes  longer  to  converge  to  the  solution  and  can  become  computationally  prohibitive. 

A  few  preliminary  tests  were  carried  out  to  pick  the  most  suitable  algorithm  for  the  multi-parameter  optimization.  The  numerical 
gradient  based  techniques  were  found  to  be  sensitive  to  function  unevenness  and  demonstrated  unstable  behavior.  The  simulated 
annealing  algorithm  was  found  robust  but  very  slow.  The  Nelder-Mead  direct  search  algorithm  provided  most  accurate  results 
for  the  MI  function  upon  reasonable  selection  of  initial  step  value  and  termination-continuation  conditions.  Briefly,  the  Nelder- 
Mead  direct  search  algorithm7  constructs  a  pyramid  in  the  six-dimensional  parameter  space  from  seven  initial  trial  points.  During 
optimization,  vertices  of  the  pyramid  converge  toward  the  solution  while  the  pyramid  shrinks  in  size.  A  new  candidate  vertex 
is  the  mirror  image  of  the  vertex  with  the  smallest  MI  value,  about  the  directly  opposite  face.  A  new  pyramid  is  built  depending 
on  the  relative  values  of  MI  at  the  candidate  vertex  and  the  other  six  vertexes.  It  is  easy  to  appreciate  that  the  optimization  will 
not  be  caught  in  a  local  maximum  with  dimensions  smaller  than  the  pyramid  size. 


A  characteristic  feature  of  ultrasound  images  is  the  presence  of  speckle  noise  that  is  likely  responsible  for  the  unevenness  in  the 
mutual  information  function.  Intra-frame  spatial  averaging,  which  also  improved  images  visually,  slightly  reduces  this  effect. 
Furthermore,  inter-frame  temporal  filtering  was  performed  that  was  particularly  effective  because  speckles  are  at  least  partially 
uncorrelated  in  adjacent  frames.  The  dynamic  nature  of  the  heart  dictates  that  only  a  few  adjacent  frames  be  averaged  however. 
Intensity  binning  when  creating  histograms  for  computing  mutual  information  was  another  method  of  speckle  reduction.  The 
net  effect  of  binning  was  to  discard  lower  voxel  intensity  bits.  The  intensity  binning  eliminated  low-power  speckles  and  their 
influence  on  MI. 


The  validation  of  our  registration  technique  requires  that  the  computed  results  be  compared  to  a  gold  standard.  The  only  gold 
standard  that  may  exist  for  the  present  application  is  the  result  of  manual  registration  by  experts.  For  the  purposes  of  this 
preliminary  study  and  to  get  an  estimate  of  the  order  of  error,  we  performed  a  set  of  intra- frameset  registration  experiemnts.  The 
orientation  of  an  ultrasound  transducer  is  virtually  fixed  during  an  examination  due  to  short  acquisition  duration.  Therefore  it 
is  appropriate  to  assume  that  the  individual  frames  of  a  frameset  are  perfectly  registered.  We  picked  two  neighboring  frames  that 
were  perfectly  aligned  to  begin  with  from  the  same  frameset,  and  applied  to  one  of  them  a  known  transformation  thus  creating 
a  known  misalignment.  The  normal  3D  spatial  registration  procedure  was  applied  to  recover  the  misalignment.  A  quantitative 
analysis  was  possible  by  comparing  the  known  transformation  with  the  solution  of  3D  spatial  registration. 


3.  DATA  PROCESSING 


3.1.  Image  Acquisition 

Ultrasound  images  used  in  the  study  were  acquired  on 
a  2.5  MHz  3D  ultrasound  scanner  manufactured  by 
Volumetries,  Inc.  The  nine  framesets  used  in  the  study 
were  images  of  the  left  ventricle  of  four  different 
patients.  The  scanner  generates  volume  frames  that  are 
sector  scans  in  3D  or  roughly  shaped  like  a  pyramid 
truncated  at  the  top  (Figure  1).  The  size  of  the  pyramid 
is  user  controlled  through  scan  depth  and  scan  angle 
parameters.  The  framesets  used  in  the  present  study  all 
had  a  scan  depth  of  16  cm  and  60  degree  of  angular  span 
in  both  B-mode  directions.  In  addition,  the  images  were 
collected  at  20  frames  per  second  and  the  number  of 
frames  in  a  set  varied  from  1 8  to  32. 


Figure  1.  The  truncated  pyramid  shape  of  3D  ultrasound 
volume  frame. 


Figure  2.  Three  orthogonal  views  showing  left  ventricle  of  a  scan  converted  3D  ultrasound  frame. 


3.2.  Scan  Conversion 

The  original  polar  data  was  scan  converted  to  a  rectilinear  grid.  The  polar  frames  had  a  64x64x512  resolution  with  512  samples 
along  the  radial  direction.  The  scan  converted  frames  were  sized  128x128x512  to  maintain  the  native  higher  radial  resolution. 
The  scan  conversion  conceptually  placed  the  smallest  bounding  box  around  the  polar  frame.  Based  on  the  value  of  acquisition 
parameters,  the  physical  dimensions  (in  mm)  of  the  bounding  box  were  167.2x  167.2x  160.  A  byte  described  the  ultrasound  voxel 
intensity  in  both  the  native  and  scan  converted  data.  Figure  2  shows  an  example  of  three  orthogonal  sections  of  a  scan  converted 
frame.  The  first  two  panels  show  two  orthogonal  long  axis  views  of  the  left  ventricle  whereas  the  third  panel  shows  a  short  axis 
view. 

3.3.  Subsampling  and  Speckle  Reduction 

The  original  128x128x512  sized  frames  were  subsampled  by  a  factor  of  two  to  keep  computational  needs  of  registration  within 
reasonable  limits.  Uniform  intra-frame  averaging  by  a  2x2x2  kernel  helped  achieve  both  subsampling  and  speckle  reduction  by 
spatial  smoothing. 


Temporal  (inter-frame)  filtering  helped  reduce  speckles  further.  A  weighted  averaging  over  three  adjacent  frames  was  performed 
with  a  [V6,  Vfe,  lA]  kernel,  i.e.,  each  voxel  was  filtered  with  analogous  voxels  from  the  preceding  and  succeeding  frames.  The 
contribution  of  each  neighboring  voxel  was  half  of  that  of  the  current  voxel.  Since  the  end  diastolic  frames  were  not  end  frames, 
the  above-described  temporal  filtering  was  carried  out  without  ambiguity.  Registration  was  performed  both  with  and  without 
temporal  filtering  to  understand  its  effects. 

An  implicit  averaging  and  hence  speckle  reduction  took  place  in  the  creation  of  histograms  for  MI  calculation.  By  binning  the 
intensities  in  the  histogram,  or  equivalently,  discarding  a  few  lower  bits  of  the  voxel  intensity,  a  further  smoothing  took  place. 
The  more  intensity  bits  were  rejected,  the  more  smoothing  took  place.  The  effect  of  intensity  binning  on  the  outcome  of 
registration  was  studied  by  repeating  experiments  six  times  by  varying  the  number  of  intensity  bits  from  all  eight  to  upper  three. 

3.4.  Data  mask 

While  the  rectilinear  representation  of  data  is  suitable  for  image  processing,  a  byproduct  of  scan  conversion  is  inclusion  of  blank 
area  outside  the  original  pyramid  into  data  representation.  The  black  region  of  zero  intensity  around  the  image  data  in  Figure  2 
corresponds  to  the  blank  area.  The  blank  area  voxels  must  be  rejected  while  computing  probability  density  functions.  To  facilitate 
this,  a  bit  mask  was  created  that  corresponded  with  the  shape  of  the  original  image  data.  This  mask  was  helpful  in  the 
determination  of  the  overlapping  zone  of  two  frames.  Although  an  analytical  formulation  for  the  region  on  overlap  is  possible, 
defining  a  mask  was  more  convenient  from  implementation  standpoint. 

3.5.  Inter-frame  MI 

As  mentioned  earlier,  adjacent  frame  MI  along  a  frameset  was  found  a  good  indicator  of  cardiac  phase.  Inter-frame  MI  was  also 
useful  in  appreciating  the  physical  meaning  of  the  mutual  information  as  a  measure  of  data  similarity.  To  this  end,  MI  was  also 
computed  between  all  possible  pair  of  frames  within  a  set. 

3.6.  Validation 

Experiments  were  carried  out  to  test  the  feasibility  of  MI  based  3D  spatial  registration  of  cardiac  ultrasound  frames  and  to 
estimate  its  accuracy  as  a  function  of  the  number  of  bits  used  for  computing  MI  and  temporal  inter-frame  filtering.  An  intra¬ 
frameset  approach  was  taken,  as  described  earlier,  where  a  known  misalignment  was  created  between  perfectly  pre-registered 
frames.  The  two  frames  selected  were  end  diastolic  frames  separated  by  four  frames  or  250  ms.  The  initial  misalignment  was 
created  by  introducing  a  shift  of  seven  units  in  each  of  the  six  transformation  parameters.  A  seven  unit  shift  in  a  translation 
parameter  was  equal  to  seven  voxels  or  approximately  9  mm.  A  seven  unit  shift  in  a  rotation  parameter  equaled  seven  degrees. 
The  sign  (positive  or  negative)  of  these  shifts  was  chosen  randomly.  Simply  reversing  the  signs  of  each  transformation  parameter 
created  a  second  initial  misalignment  between  a  pair  of  frames.  Overall,  3D  spatial  registration  was  performed  on  nine  frame 
pairs  originating  from  nine  framesets  and  each  frame  pair  was  registered  twice  from  two  different  starting  misalignments. 
Furthermore,  registration  was  launched  six  times  for  each  starting  misalignment  to  study  the  relationship  of  the  accuracy  of 
registration  with  the  number  of  intensity  bits  used  for  computing  MI. 

If  the  registration  worked  perfectly,  a  pair  of  frames  would  have  come  back  to  the  original  alignment  upon  registration.  In  reality, 
the  algorithm  converged  to  a  solution  that  did  not  fully  annul  the  user-introduced  misalignment.  The  difference  in  the  expected 
and  the  actual  solutions  was  our  evaluation  criterion.  This  difference  was  represented  as  a  gross  displacement  and  angle 
mismatch. 


4.  RESULTS 

The  adjacent  frame  MI  versus  the  acquired  ECG  for  two  different  frame  sets  is  shown  in  the  two  panels  of  Figure  3.  The  midpoint 
of  end  diastole  is  typically  about  200  ms  left  of  the  peak  of  R  wave.  The  maximum  of  adjacent  frame  MI  approximately  coincided 
with  end  diastole.  The  low  values  along  the  MI  curve  depict  phases  where  the  movement  of  the  heart  between  frames  was  more 
drastic.  Another  way  to  present  the  same  concept  is  a  surface  plot  of  inter-frame  MI  between  all  pairs  of  frames  in  a  frameset 
(see  Figure  4).  The  frameset  used  for  this  illustration  contained  32  frames.  The  tallest  ridge  corresponds  to  MI  of  a  frame  with 


itself.  The  other  peaks  correspond  to  pair  of  frames  with  end  diastolic  phase.  It  should  be  noted  that  the  frameset  used  to  create 
the  plot  here  spanned  roughly  two  cardiac  cycles  and  hence  two  end  diastolic  phase  frames  explaining  the  presence  of  multiple 
peaks. 


Figure  3.  Two  examples  of  "adjacent  frame"  mutual  information  versus  ECG. 


Figure  4.  Mutual  information  between  all  pairs  of  32  frames  of  a  frameset:  The  diagnonal  ridge  corresponds  to  MI  of  a 
frame  with  itself.  Other  peaks  correspond  to  pair  of  frames  with  end  diastolic  phase. 


Table  1  compiles  the  result  of  registration  for  three  randomly  selected  pairs  of  frames.  The  displacement  and  angle  mismatches 
from  the  expected  solution  is  presented  for  six  varying  schemes  of  MI  computation.  The  optimal  number  of  intesity  bits  to  use 
for  MI  computation  is  not  evident  from  the  three  cases  presented.  A  more  comprehensive  trend  emerges  from  the  complete  results 
enumerated  below. 


Table  1.  The  accuracy  of  registration  for  three  pairs  of  frames  for  six  different  schemes  of  computing  MI. 


Frameset  ID  and 
frame  serial  numbers 

Number  of  intensity  bits 

8 

7 

6 

5 

4 

3 

BaaEH 

— 

1 

[D802, 

frames  17  &  12] 

0.36  mm 
1.19  deg 

— 

— 

HHi 

KE9RH 

— 

ME3EHH1 

■m 

— ■ 

1.23  mm 
1.00  deg 

0.81  mm 
0.09  deg 

0.89  mm 
0.38  deg 

Results  of  registration  for  all  18  combinations  of  nine  framesets  and  two  initial  misalignments  are  presented  in  Table  2.  Each 
row  has  results  specific  to  the  number  of  intensity  bits  used  for  MI  calculation.  Both  root-mean-square  (RMS)  and  median 
displacement  and  angle  values  are  reported.  The  median  values  are  reflective  of  the  most  common  case  without  the  influence 
of  outliers  that  the  RMS  values  suffer  from.  The  last  column  has  the  number  of  outliers,  which  were  the  cases  for  which  the  result 
was  off  of  the  expected  solution  by  more  than  the  cutoff  values  of  6  mm  and  6  degrees.  The  quantity  "Deviation  from  maximal 
MI"  measures  the  distance  of  search  results  for  a  given  number  of  bits  from  the  result  associated  with  the  largest  MI  upon 
convergence.  This  quantity  helps  us  appreciate  the  spread  of  the  search  results.  The  RMS  and  medial  values  of  displacement  and 
angle  with  respect  to  maximal  MI  solution  are  tabulated  in  third  and  fourth  columns  of  Table  2. 


Table  2.  Average  mismatch  upon  convergence  for  six  schemes  of  computing  MI. 


Number  of 
intensity 
bits 

RMS  deviation  from 
perfect  solution 
(mm;  degree) 

Median  deviation 
from  perfect 
solution 
(mm;  degree) 

RMS  deviation  from 
maximal  MI 
(mm;  degree) 

Median  deviation 
from  maximal  MI 
(mm;  degree) 

Number  of 
outliers 

8 

3.62;  6.98 

1.40;  2.54 

3.36;  6.88 

0.58;  2.11 

7 

7 

1.84;  4.02 

1.27;  1.02 

1.59;  3.93 

0.81;  0.96 

5 

6 

1.84;  1.04 

1.08;  0.75 

1.45;  0.95 

0.56;  0.68 

0 

5 

1.87;  2.45 

1.27;  0.95 

1.46;  2.40 

0.51;  0.86 

1 

4 

2.01;  2.79 

1.43;  0.59 

1.60;  2.70 

0.80;  0.77 

2 

3 

2.23;  4.22 

1.16;  0.77 

1.95;  4.19 

0.54;  0.71 

3 

Table  3.  Average  mismatch  upon  convergence  for  six  schemes  of  computing  MI  for  temporally  filtered  frames. 


Number  of 
intensity 
bits 

RMS  deviation  from 
perfect  solution 
(mm;  degree) 

Median  deviation 
from  perfect 
solution 
(mm;  degree) 

RMS  deviation  from 
maximal  MI 
(mm;  degree) 

Median  deviation 
from  maximal  MI 
(mm;  degree) 

Number  of 
outliers 

8 

2.17;  4.73 

1.47;  0.70 

1.85;  4.71 

0.35;  0.59 

3 

7 

1.67;  2.40 

1.16;  0.51 

1.27;  2.38 

0.24;  0.45 

2 

6 

1.57;  0.86 

1.02;  0.70 

1.18;  0.77 

0.30;  0.29 

0 

5 

1.64;  0.78 

1.03;  0.51 

1.26;  0.72 

0.35;  0.55 

0 

4 

1.63;  0.85 

1.09;  0.54 

1.27;  0.80 

0.29;  0.41 

0 

3 

1.77;  2.41 

1.18;  0.77 

1.37;  2.38 

0.51;  0.68 

1 

The  same  results  as  presented  in  Table  2  are  presented  in  Table  3  with  the  difference  that  registration  was  performed  on 
temporally-filtered  frames.  Improvement  in  the  accurary  of  registration  is  obvious.  The  number  of  outliers  also  went  down 
significantly.  From  Table  2  and  Table  3,  is  evident  that  using  either  5  or  6  intensity  bits  to  compute  MI,  together  with  temporal 
filtering,  produces  most  accurate  results. 

It  deserves  to  mention  that  the  value  of  MI  for  registered  frames  for  all  nine  pairs  ranged  between  0.3  and  0.95  when  using  all 
8  bits,  and  0.1 1  and  0.30  when  using  3  bits.  This  shows  data  dependency  of  MI  function.  An  unexpected  observation  was  a 
slightly  larger  MI  value  at  the  point  of  convergence  than  that  at  the  zero  point  or  perfect  solution  almost  all  the  time.  This 
demonstates  that  the  search  algorithm  behaved  satisfactorily  and  that  the  unevenness  of  the  MI  function  prevented  the  search 
to  converge  to  the  perfect  solution. 

Execution  time  for  each  registration  equaled  the  product  of  time  for  the  evaluation  of  MI  and  the  total  number  of  these 
evaluations  along  the  search  path.  On  an  average,  evaluation  of  MI  took  approximately  8  seconds  on  a  450  MHz  Pentium-II 
processor  computer.  Most  of  this  time  was  spent  in  trilinear  interpolation  of  the  secondary.  Typically  there  were  150-250  such 
evaluations  and  therefore  the  total  processing  time  was  between  20-30  minutes. 

Inter-frameset  registration,  an  example  of  which  is  shown  Figure  5,  was  performed  although  not  quantitatively  validated.  The 
figure  shows  three  pairs  of  long-axis  views.  The  middle  pair  shows  two  views  of  the  primary.  The  left  pair  shows  the  identical 
views  of  the  secondary  before  registration.  The  right  pair  shows  the  same  views  of  the  secondary  following  registration.  A  visual 
inspection,  at  least  qualitatively,  demonstrates  the  success  of  inter-frameset  registration. 


(a)  (b)  (c) 

Figure  5.  Example  of  inter-frameset  registration:  (a)  secondary  frame  before  registration,  (b)  fixed  primary  frame,  and  (c) 
secondary  frame  after  registration.  The  x-z  plane  cross-sections  are  in  the  upper  row  and  the  y-z  plane  cross-sections  in 

lower. 


5.  DISCUSSION  AND  CONCLUSIONS 


We  have  demonstrated  that  mutual  information  based  registration  of  ultrasound  images  is  feasible.  We  have  also  presented  a 
procedure  to  compute  mutual  information  such  that  the  success  and  accuracy  of  registration  is  maximized.  A  quantitative 
validation  of  registration  was  performed  using  two  end  diastolic  frames  from  the  same  frameset.  The  overall  accuracy  was  of 
the  order  of  a  voxel  dimension.  The  present  study  provides  a  good  foundation  for  further  investigation  and  quantitative  validation 
of  inter-frameset  registration  and  its  application  to  real  life  applications. 

The  computation  of  mutual  information  is  dependent  on  individual  and  joint  histograms  that  in  turn  depend  on  voxel  intensities. 
Presence  of  speckles  in  ultrasound  images  appears  to  cause  unevenness  in  the  mutual  information  function.  This  unevenness 
poses  two  problems:  firstly,  it  can  hide  or  diminish  the  desired  global  maximum  and,  secondly,  it  creates  several  local  maxima 
that  may  confuse  the  search  algorithm.  An  evidence  of  this  was  seen  when  the  MI  function  did  not  have  a  global  maximum  at 
perfect  solution.  Therefore  the  search  algorithm  functioned  as  expected  and  the  error  can  be  attributed  to  MI  function. 

Temporal  filtering  was  effective  in  creating  a  smoothly  varying  MI  hypersurface  as  the  accuracy  of  registration  improved 
significantly  with  it.  The  intensity  binning  also  had  a  significant  effect  on  the  shape  of  MI  hypersurface  and  the  overall  accuracy. 
The  experiments  suggested  that  using  five  or  six  upper  bits  of  intensity  to  compute  MI  was  most  effective.  A  possible  explanation 
is  that  using  all  eight  bits  preserves  the  speckle  noise  what  corrupts  the  underlying  signal.  As  the  lower  bits  were  discarded 
gradually  signal-to-noise  ratio  (SNR)  first  increased  due  to  reduction  in  noise  and  then  started  to  decrease  due  to  faster  loss  of 
signal  than  noise.  The  SNR  peaked  for  five  or  six  bits  of  intensity.  Stated  differently,  at  intermediate  levels  (5  or  6  bits)  the  MI 
function  is  based  on  the  most  salient  features  of  the  histograms.  At  higher  levels,  noise  starts  interfering  with  the  features  whereas 
at  lower  levels  salient  features  start  to  wash  away. 

The  mutual  information  based  registration  does  not  have  an  infinite  capture  range  (i.e.,  the  maximum  misalignment  from  which 
the  search  algorithm  can  still  converge  to  the  desired  solution).  If  the  initial  misalignment  is  large,  the  convergence  may  not  take 
place.  A  misalignment  of  seven  units  was  large  enough  to  be  observable  and  yet  small  enough  to  be  within  the  capture  range  for 
our  datasets.  Since  volume  of  overlap  factors  into  the  MI  calculation,  the  MI  function  can  have  very  high  values  for  very  small 
overlap.  Quite  clearly,  those  values  do  not  correspond  to  the  desired  solution  where  a  considerable  overlap  is  expected.  The 
search  algorithm  seeks  to  converge  to  the  global  maximum  within  this  capture  range.  The  capture  range  is  also  dependent  on  the 
size  of  the  initial  pyramid  in  Nelder-Mead  algorithm.  However,  the  size  cannot  be  made  arbitrarily  large  without  making  the 
algorithm  unstable. 

The  intra-frameset  validation  was  a  good  indicator  of  the  accuracy  of  the  overall  process.  When  MI  function  was  computed  in 
the  most  optimal  way,  the  typical  error  was  of  the  order  of  a  voxel  dimension.  Achievement  of  subvoxel  accuracy  was  beyond 
the  scope  of  this  study.  Although  we  have  achieved  visually  acceptable  results  of  inter-frameset  registration,  a  quantitative 
validation  with  expert  data  is  desired. 

Encouraged  by  the  results  of  this  preliminary  investigation,  our  future  work  will  focus  on  inter-frameset  registration.  We  will 
investigate  ways  to  compute  MI  function  for  further  accuracy,  perform  rigorous  validation,  lower  execution  time  and  attempt 
clinical  deployment  of  the  developed  methods. 
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